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ABSTEACT 


This  study  concerns  the  laminar  flow  of  an  incompressible  fluid  through 
a  curved  pipe  with  an  elliptical  cross  section.  The  governing  equations  are 
derived  by  applying  the  Navier-Stokes  and  continuity  equations  in  cylindrical 
coordinates  and  the  method  of  successive  approximations  to  get  the  five  partial 
differential  equations.  These  equations  are  solved  by  the  perturbation  method, 
and  twenty  numerical  examples  are  presented.  For  each  example,  arbitrary 
numerical  parameters  are  assumed  for  input  into  an  IBM  7094  Computer  to 
obtain  solutions  for  the  simultaneous  algebraic  equations.  For  an  eccentricity 
of  one,  the  ellipse  degenerates  to  a  circle  and  Dean's  solution  for  the  stream¬ 
line  flow  of  an  incompressible  fluid  through  a  curved  pipe  with  a  circular 
cross  section  is  obtained. 
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CHAPTER  I 

INTRODUCTION 

The  design  of  present  and  future  high  performance  missiles,  rockets, 
and  aircraft  necessitates  the  optimum  use  of  all  available  space.  To  help 
accomplish  u..s  requirement,  bent  tubing  and  curved  pipes  are  frequently  used. 
In  the  design  of  these  piping  systems,  it  is  gen  erally  assumed  that  the  circular 
cross-sectional  area  remains  constant  to  simplify  the  analysis.  This  is  a  true 
assumption  only  if  the  pipe  is  cast  or  forged,  or  if  a  fitting  with  a  circular 
cross  section  is  used  to  obtain  the  required  curvature.  However,  a  fitting  is 
used  only  for  small  curvatures.  Large  curvatures  are  usually  desired  to 
reduce  the  pressure  drop,  decrease  the  size  of  pipe  or  increase  the  flow  rate. 
This  is  accomplished  by  bending  a  piece  of  straight  pipe  to  the  desired  curva¬ 
ture,  thus  distorting  the  circular  cross  section  and  making  it  elliptical. 

The  purpose  of  this  report  is  to  study  the  effect  of  the  eccentricity  of 
an  elliptical  cross  section  in  a  pipe  with  a  small  curvature  and  to  get  a  second- 
order  approximation  of  the  flow  rate  through  the  pipe. 

The  term  "curved  pipe"  as  used  in  this  report  is  a  pipe  with  an  elliptical 
cross  section  bent  so  that  the  center  of  the  ellipse  forms  a  circular  arc.  It  is 
assumed  that  the  flow  is  fully  developed  throughout  the  region  and  that  the 
straight  pipes  are  connected  to  the  bend  so  that  they  always  lie  in  the  plane  of 
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the  bend  with  their  center  lines  tangent  to  the  center  line  of  the  bend. 

Dean  [1927]  showed  that  the  motion  of  a  viscous,  incompressible  fluid 
in  a  curved  pipe  with  a  circular  cross  section  consists  of  a  primary  motion 
along  and  parallel  to  the  center  line  of  the  pipe,  and  the  secondary  motion 
which  is  in  the  plane  of  the  cross  section.  His  solution  is  the  most  detailed 
and  perhaps  the  best  theoretical  study  performed  to  this  date  although  it  is 
severely  limited  to  Dean's  numbers  less  than  about  400.  This  results  from 
neglect  of  the  unsymmetrical  terms  in  the  second-order  approximation  of  the 
flow  rate  and  the  insufficient  number  of  terms  used  with  the  method  of  suc¬ 
cessive  approximations.  He  also  assumed  small  pipe  curvature  and  that  the 
circular  cross  section  remained  circular  after  the  pipe  was  bent.  Within  its 
limitations,  Dean's  solution  agrees  very  well  with  experimental  data,  and  his 
theory  has  been  widely  used  and  extended. 

Thomas  and  Walters  [1963]  used  Dean's  approach  and  analyzed  his 
problem  using  an  elastico-viscous  liquid.  Their  work  showed  that  the  e1ostico- 
viscous  property  of  the  fluid  reduced  the  curvature  of  the  streamlines  in  the 
central  plane  and  increased  the  rate  of  flow  through  the  pipe. 

Baura  [1963]  took  Dean's  [1928]  basic  equations  and  integrated  across 
the  boundary  layer  to  get  the  momentum  integral  form.  The  equations  were 
than  solved  by  Polvhausen's  method.  The  results  agreed  very  well  with 
experimental  data  but  are  probably  no  more  accurate  than  Dean's  solution. 
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Clegg  and  Power  [1963]  analyzed  the  flow  of  a  Bingham  fluid*  in  a 
slightly  curved  pipe,  but  their  solution  was  accurate  only  to  the  first-order 
approximation  and  did  not  include  the  effect  of  the  curvature.  The  effects  of 
a  plug  being  inserted  in  the  center  of  the  pipe  were  also  studied. 

Thomas  and  Walters  [1965]  studied  the  flow  of  a  fluid  through  a  curved 
pipe  with  an  elliptical  cross  section,  but  their  results  were  accurate  only  to  the 
first-order  approximation.  Their  solution  indicated  that  the  rate  of  flow 
through  the  pipe  was  independent  of  the  curvature.  It  was  then  concluded  that 
the  second-order  terms  were  required  to  determine  the  effect  of  curvature  on 
the  rate  of  flow  through  the  pipe. 

The  experimental  approach  to  special  cases  of  this  problem  was  taken  by 
White  [1929]  and  Keuleganand  Beij  [1937].  The  results  obtained  conformed  very 
well  with  Dean's  solution.  White  performed  experiments  with  water  and  oil 
flowing  through  curved  pipes  with  oval  and  circular  cross  sections  to  determine 
the  law  of  resistance  for  streamline  flow.  Keulegan  and  Beij  conducted  experi¬ 
ments  with  water  flowing  through  a  curved  pipe,  with  a  circular  cross  section 
prior  to  bending,  to  determine  the  pressure  losses  ^ue  to  the  large  curvature  of 
the  pipe.  The  first  experiments  resulted  in  an  eq  Jon  for  the  prediction  of 
head  loss,  and  the  second  experiments  resulted  in  an  equation  for  the  increase 
in  resistance  in  a  bend  as  compared  to  that  of  a  straight  pipe. 

This  study  considers  the  streamline  motion  of  a  viscous,  incompressible 

*  A  Bingham  fluid  is  a  material  that  can  support  a  finite  stress  elasti¬ 
cally  without  flow  and  flows  with  constant  plastic  fluidity  when  the  stresses  are 
sufficiently  great. 
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fluid  in  a  curved  pipe  with  an  elliptical  cross  section.  The  angle  of  bend  in  the 
pipe  is  unrestricted  provided  that  the  curvature  is  small  and  the  flow  in  the 
bond  is  fully  developed. 

The  governing  partial  differential  equations  are  derived  in  Chapter  II 
by  applying  the  equations  of  motion  ( Navier-Stokes  equations)  and  the  continuity 
equation  in  cylindrical  coordinates.  The  method  of  successive  approximations 
is  then  used  to  obtain  equations  which  are  accurate  to  the  second-order  approxi¬ 
mation. 

In  Chapter  III  these  nonlinear  partial  differential  equations  are  solved. 
Simultaneous  equations  involving  the  constants  are  then  obtained  by  matching 
coefficients.  The  unsymmetrical  terms,  although  negligible  for  small  Dean’s 
numbers,  are  not  small  for  large  Dean’s  numbers  and  are  not  neglected  in  the 
solutions.  Therefore,  these  solutions  are  more  accurate  than  the  corresponding 
second-order  approximation  solution  presented  by  Dean  for  a  pipe  of  circular 
cross  section.  By  letting  the  cross-sectional  ellipse  degenerate  to  a  circle,  the 
equations  of  the  present  work  reduce  to  Dean's  [1928]  equations  except  for  the 
last  equation,  in  which  Dean  neglected  the  unsymmetrical  terms. 

The  four  sets  of  simultaneous  equations  in  matrix  form  are  presented  in 
Chapter  IV,  and  their  computer  solutions  are  in  Appendix  B.  The  equations  for 
the  rate  of  flow  through  the  curved  pipe  are  integrated,  and  equations  for  20 
different  eccentricities  of  the  cross-sectional  ellipse  are  presented  in  Chapter  FV. 

The  discussion,  summary  and  conclusions,  and  recommendations  for 
future  research  are  given  in  Chapters  V,  VI  and  VII,  respectively. 
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CHAPTER  n 

GOVERNING  DIFFERENTIAL  EQUATIONS 

A.  Assumptions 

This  work  is  a  stud}7  of  the  streamline  motion  of  a  Newtonian  fluid  in  ?. 
curved  pipe  with  an  elliptical  cross  section.  The  equations  of  motion  in  cylin¬ 
drical  coordinates  (r\  9,  y ’ )  are  applied  to  the  arrangement  shown  in 
Figure  1. 

The  following  assumptions  are  made  in  the  derivation  of  the  governing 
partial  differential  equations : 

1.  Streamline  motion  theory  applies. 

2.  Body  forces  are  negligible. 

3.  The  flow  is  steady,  uniform  and  incompressible. 

4.  The  curvature  of  the  pipe  is  small. 

5.  The  flow  is  axisvmmetric  and  all  flow  variables  except  pressure 
are  independent  of  0. 

6.  The  flow  is  fully  developed  throughout  the  region  of  study. 

For  laminar,  incompressible  flow  the  unknown  variables  consist  of  the 
pressure  and  three  components  of  velocity.  To  solve  for  these  four  unknowns 
the  three  equations  of  motion  (Navier-stokes  equations)  and  the  continuity 
equation  are  used.  These  equations  appear  in  Schlichting  |1968]  as  well  as 


Figure  1.  Cylindrical  Coordinate  System 
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most  text  books  on  fluid  mechanics,  and  are  given  in  the  following  sections. 

B.  Navier-Stokes  Equations  in  Cylindrical  Coordinates 
For  the  r' -direction: 


aqr  +  ,%d_%+  d_\ 

at  +  qr  9r'  T  r'  a©  ay'  “  r' 


(2. 1) 


ap  (  2  qr  2  9q© 

r  8r’  \  *  — r2  „f2  80 


For  the  0  -direction: 


9qe  aq©  q09q© 

p|__-  +  q  — —  +  q  —  + 


at 


?r  ar'  +  r'  a©  "  ‘V  9y* 


8q©  qr<fe\ 


n 


f  _  1  ap  .  Lyl  2  8q*  q° 

pf0  r’  a©  <3*e+  ^,"3© 


(2.  2) 


For  the  y' -direction: 


fdq  3q  q  aq_  do 
y  y  ^0  v  y 

p  — +  q  — ^  + —  — -  +  q  — £- 

\  9t  4r  3r'  r'  30  y  dy' 


(2.3) 


ap 

y  9y' 


=  Pf~  ~  ^  + 


C.  Continuity  Equation  in  Cylindrical  Coordinates 


3q  q 
r  r 

■■■  •+• 


:  %  +  '% 

r'  '  r'  90  3y' 


(2.4) 
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where  f  ,  f_ ,  and  f  are  the  components  of  the  body  force;  q  ,  q  ,  and  q  are 
r  0  y  r  o  y 

the  components  of  velocity;  P  is  the  pressure;  p  is  the  fluid  density;  p  is  the 


a2  1  9  i  a2  a2 


absolute  viscosity  of  the  fluid  and  V2  =  — +  —  — —  + — j  -rj  + - -  . 

r’  9r'  r'  90*  „  ,2 

ar'2  9y 


The  assumption  of  steady  flew  and  negligible  body  forces  reduces 
equations  (2.1)  through  (2.3)  to: 


( 

“  Vr  ar’  +  r'  30 


S'  2  8qel 


+  jp  +  ,v,!q_ 

y  9y'  r-  ar'  ^  .2  .2  9© 


(2.5) 


r'  r' 
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9q  q„  9q  ?c  , 
y  0  y  q  y\ 

p  in  — ■*- + - *-  +  v — *■] 

'4r  3r’  r'  80  y  9y' 


_!£+  Mv’2q 
9yf  y 


(2.7) 


Since  the  flow  is  uniform  and  q  ,  q  ,  and  q  are  independent  of  0, 

r  y  v 


aqr  9qy  9q0  p 

— —  =  — ^-=  — —  =  0.  The  kinematic  viscosity,  v  -  ~  and  the  equations 
80  90  80  ”  P 


<2.  4)  through  (2.7)  become: 


8q  q  8q 
r  r  v 
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8r'  r'  9y' 


[2.8) 
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Figure  2.  Cartesian  Coordinate  System  for  Any  0 


(2. 12) ,  (2. 13) ,  and  (2. 14)  are  fuictions  ofr'  and  y'  only  except  the  pressure 

i  a  (P/p)  d(P/p)  ,  3(P /p)  _  ,  .  . 

terms  -  —  - , - - — ,  and  -  — — —  .  Solving  each  of  these  equations 

R  3©  3r  3y 

for  its  pressure  term  shows  that  the  three  pressure  terms  are  also  functions  of 


r'andy1  only.  Therefore,  — =  e,©  +  f(r\y’),  where  et  is  a  constant.  Then 


>  .,(7) 

R  3© 


1  iL 

R  3© 


fe,©  +  f(r’,y’ •  ] 


=  -  — =  constant  .  (2,15) 


0  G 

Redefining  the  constant  — —  =  — makes  G  a  constant  which  can  be 

n  p 

called  the  mean  pressure  gradient.  It  is  the  space-rate  of  decrease  in  the 
pressure  along  the  central  line  traced  out  by  the  center  of  the  pipe. 

Differentiating  equations  (2.12)  and  (2.14)  with  respect  toy’  andr’, 
respectively,  and  subtracting  (2.14)  from  (2.12)  eliminates  the  pressure 
terms.  The  results  are 


3q 
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3r 


r  3  r 


r  -  q. 


\  + 

v  By'  /  3y’  \^r  3r'  3v' 
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(2.16) 
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Using  equation  (2.15)  and  -  =  ■—  .equation  (2.13)  can  be  expressed 


as 
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3q, 


3r’ 


y  9y’ 


©  + 
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~+  "(v-V  -  ‘e 


0 


(2.17) 


Because  of  the  highly  nonlinear  character  of  these  equations,  a  solution 
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cannot  be  obtained  in  this  form.  Tt  is,  therefore,  necessary  to  simplify  them 
in  such  a  way  that  the  nonlinear  and  curvature  effects  are  retained  To  accom¬ 
plish  this  and  transfer  the  origin  from  O'  to  0  (Figure  1)  it  is  assumed  that  V’2 

2  8  9  1 

in  cylindrical  coordinates  equals  V'2  in  Cartesian  coordinates,——  =  — -  + — — 

3x  9r  r 

and  r'  is  approximately  R.  The  small  curvature  of  the  pipe  provides  the  basis 
for  these  assumptions.  It  is  evident  that  these  assumptions  are  equivalent  to 
the  ones  made  by  Dean  [1928]  for  a  curved  pipe  with  a  circular  cross  section 
because  the  equations  obtained  degenerate  to  those  presented  by  Dean  [1928]. 
for  the  circular  case.*  Since  the  i '-  and  x'-directions  are  the  same,  q^  =  q^, 
where  q^  is  the  velocity  in  the  x'-direction.  Use  of  tnese  assumptions  reduces 
equations  (2.8),  (2.16),  and  (2.17)  to: 


q  9q 

q  +  _2L  +_X,  o 

[ax'  r  r'  9y' 


(2.13) 
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where,  by  definition,  V“ 


9x 


,2  +  ,2  in  the  Cartesian  coordinate  system. 

ay 


*The  only  exception  is  Dean's  solution  for  the  second-order  approxima¬ 
tion  of  the  velocity  in  the  9-direction,  in  which  he  neglected  the  unsymmetrical 
terms  of  the  solution. 
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U  = 


d<P 

dy 


V  =  - 


alp 

ax' 


(2.33) 


Substituting  equations  (2.33)  into  equation  (2.29)  yields: 


dx  8y  dy  d x  r 


d2tp  dip  dip  d2ip 
^2  dy  3x  dxdy 


+  KW 


aw 

ay 


V2(vV) 

r  ay2 


(2.34) 


The  terms  — 
r 


d2ip  dp  dp  d2ip 


ax' 


2  ay  ax  dxdy 


i  afy 

and - are  third-order  terms  and  are 

r  ay2 


negligible  for  a  solution  accurate  to  the  second  approximation.  Therefore, 
equation  (2.34)  is  expressed  as: 

/3*  a  aif  ^  tsw™j£  +  i! V*  .  (2,35) 


ax  ay  ay  ax 


ay 


,ax2  av2y 


Substituting  equations  (2.33)  into  equation  (2.30)  yields: 


it  aw  _  aw  =  c  +  v,w 

ay  ax  ax  ay 


(2.J6) 


The  asymptotic  expansions  used  to  apply  the  method  of  successive 
approximations  are  as  follows: 

ip  =  Kipt  +  K2ip2  +  ...  (2.37) 


W  =  W0  +  KW,  +  K2W?  +  .  . . 


Different,  uting  equations  (2.37)  and  substituting  into  equation  (2.35)  yields: 


16 


K  *>  +  K*  «A»  -  (k  *!  +  Ks  5&L 


8x 


3x  |av 


9y 


ay  /  ax 


[k 


(2.37, 


ax* 


ax2 


+  K^+K^+. 


ay' 


K(W0 


ar 


+  KW,  +  K2W2)^li 


+  x 


»«!  ♦**.».  -(£L+.jLU*u  k* 


vy 


ay 


yax2  ay2/  \  axz 


K»V,+  k-2  ^*2 


av 


.2 


o 

av^ 


9^2 


>x 


Equating  the  zero-order  terms  (involving  K°): 


0=0 


Equating  the  first-order  terms  (involving  K1) : 


vVi  =  ~w0 


aw0 
ay  ’ 


38) 


where  V4#t  =  by  definition  and 


v2  =  Z_  +  .al 

Sx2  3y2 

Equating  the  second-order  terms  (involving  K2): 


X7% 


dtp-i  d3il> i  a^j  a3''i  a^'i  aV3  a</q  aty, 
^  IT"  ^  axTy2  "  75  a"7Ty  "  ^  ^ 


(2.39) 


aw, 


awn 


-w«lf "  WllF 


Differentiating  equations  (2.37)  and  substituting  into  equation  (2.38)  yields: 


> 
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/  \ 


-  K 


+  k2  \  +  K  +  y}  aM 

ax  ax  ]\  3y  9y  Sy  / 


_  /  9?w0  „  anv,  2  92w2  a  w, 

=  C+( -  +  K - L+Kz».  i+  - £ 


\  axz 


ax2 


+  k!^  +  k’3.5M  . 


3y 


ay 


ay 


Equating  the  zero-ox der  terms  (inv -’ng K°) : 


V2W0  =  -c 


(2.40) 


Equating  the  first-order  terms  (involving  K1) : 


v2w  =  9Wo  _  aw0 

1  9y  ax  ax  ay 

Equating  the  second-order  terms  (involving  K2): 


(2.41) 


V2w  =  +  ail  9Zi  -  1  dli-dJldjy±.  (2  40) 

2  ay  ax  'ay  ax  ~ax  ay  ax  ay  * 

Terms  involving  K  to  powers  higher  than  two  are  negligible  because  the  equa¬ 
tions  are  only  accuJ.rte  to  the  second-order  approximation. 

Equations  (2.38)  through  (2.42)  are  in  a  solvable  form  and  are  the 
governing  partial  differential  equations.  The  boundary  conditions  for  equations 
(2.38)  and  (2. 39)  are: 


=  ij>2  =  0  when  1  -  x2  -  mV  -  0 


(2.43) 
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and 


dyt  _  8e2  _  a^j  _  d4> 2  _  Q 

ax  "  ax  “  “ay-  W  _ 


■when  1  -  x2  -  mV2  =  0 


(2.43) 


where  m  is  a  positive  constant  related  to  the-  eccentricity  of  the  ellipse,  e,  by 
the  equations 


The  boundary  conditions  for  equations  (2.40)  through  (2.42)  are: 

W0  =  W,  =  W2  =  0  ,  when  1  -  x2  -  m2y2  =0  .  (2.44) 

2 

Equation  (2.40) ,  V  \V0  =  -C,  is  the  governing  partial  differential 
equation  for  flow  through  a  straight  pipe.  This  is  the  same  equation  presented 
by  Dean  [1928].  Therefore,  the  boundary  value  nroblem  for  flow  through  a 
straight  pipe  with  a  circular  cross  section  differs  with  that  for  a  straight  pipe 
with  an  elliptical  cross  section  only  in  the  boundary  conditions. 
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CHAPTER  ID 

SOLUTIONS  TO  EQUATIONS 

The  nature  of  differential  equations  (2,38)  through  (2.42)  and  the 
geometry  of  the  cross  section  of  the  pipe  suggest  that  the  solutions  for  W0,  Wj, 
\Y;,  vj.  and  i2  may  be  assumed  to  be polynomials  in  x  and  y.  For  convenience, 
the  solutions  are  assumed  in  a  form  so  that  the  symmetry  requirements  and 
boundary  conditions  are  satisfied  in  advance.  The  boundary  conditions  on  d} 
and  i2  are  given  by  equation  (2.43) ,  and  the  condition  on  W0,  W,f  and  W2  is  that 
each  vanishes  on  the  boundary  as  required  by  equation  (2.44) .  Inspection  of 
the  governing  differential  equations  shows  that  W0  and  W2  are  symmetric  in 
both  x  and  y;  Wi  is  symmetric  in  y  but  antisymmetric  with  resoect  to  x:  is 
symmetric  in  x  but  antisymmetric  in  y:  and  C2  is  antisymmetric  with  respect  to 
both  x  and  y.  Assumed  so. ...ions  which  satisfy  these  requirements  are: 


W0  =  All  -  x2  -  m'V)  (3.1) 

C]  =  (1  -  X2  -  m2y2)y(A0  +  AjX2  +  A,y2)  (3.2) 

Wj  -  (1  -  x2  -  m2y2)  x  (b0  ~  b,x2  b2y2  +  fyx4  +  fyx2y2  +  bgy4 

(3.3) 

-1-  bgx6  +  b7x4y2  +  bgx2y!  -  bHye) 

d 2  =  (1  -  x2  -  m2y2)2x  y  (C0  +  C}x2  +  C2y2  +  C3x4  +  CjX2y2 

(3-4) 

+  Cc,y4  +  Cgx6  +  C7xsy2  +  C8x2y4  -  C,/  +  C10x8 


+  C,jxV  4  C12x4y4  -  Cnx2yfi  +  Cuy8) 
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W2  =  (1  -  x2  -  mV)  (dc  +  djX2  *  djy2  +  d3x4  +  d4x2y2  +  d^4  (3.5) 

+dsxe  -  d7xy2  +  dgx2y4  -  d9y6  +  d]{>xa  +  d^xV  +  d12x4y4 

-  di.ixV  *  di4y8  +  disx'°  +  dt6xV  -  dnx8yJ  +  di8x4y6 

-  d19x2y8  +  d20y!0  -K  d21x12  +  d22x10y2  +  d2,x8y4  -  d24xeys 
-d25xJy3  -  d26x2y,°  *  d2Tv12  +  d28x14  f  d,,x12y2  +  d30x10y4 
*  dnA6  +  d:i2x6y8  -  djjxV0  -  d31x2y12  +  d35yn) 


The  coefficients  of  the  polynomial  terms  in  these  solutions  will  be  found  by 


substituting  them  into  the  governing  differential  equations  and  by  matching  the 
coefficients. 

Partially  differentiating  equation  (3. 1) .  substituting  into  equation  (2.40) 
and  sohing  for  A  gives 


Partial  differentiation  of  equations  (3.1)  through  (3.5),  substituting 
the  expressions  into  each  of  equations  (2.38)  through  (2.42)  and  equating  the 
coefficients  of  like  terms  reduces  the  problem  to  the  solution  of  four  sets  of 
simultaneo  s  equations,  each  involving  the  constants  in  the  corresponding- 
assumed  solution. 

Partially  differentiating  equations  (3.1)  and  (3.2),  substituting  into 
equation  (2.38)  and  equating  the  coefficients  of  like  terms  yields: 


) 
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a2  2 

(5m4  +  2m2  +  1)A0  +  (-2m2  -  2) A,  +  (-10m2  -  2)A;  =  — (3.7) 

a2-,4 

(10m4  +  6m2)  A,  +  (105m4  +  20n2  +  3)A2  =  -—f -  (3.8) 

A2  2 

(5m4  +  12m2  +  15)  A,  +  (10m2  +  6)A2  - - -  (3.9) 

Partially  differentiating  equations  (3.1)  through  (3.3),  substituting  into 
equation  (2.41)  and  equating  the  coefficients  of  like  terms  yields: 

'm2  -  3)bo  -  3b]  -  b2  =  AA«  (3.10) 

(m2  +  10)  bj  -  bo  -  10  b3  -  b4  =  AAj  -  2AA0  (3.11) 

3m2bj  -  (6m2  -  3)b2  -  3b4  -  6bs  =  3AA2  +  4AA0m2  -  6AA0m2  -  2AAtm2 

(3.12) 

(m2  -  21) b3  ~  b4  -  21b6  -,b7  -  AA„  -  2AA]  (3.13) 

5m2b3  -  (3m2  -  o)b4  -  3b;  -  5b7  -  3b8  =  AA0m2  -  AA4m2  -  3AA2 

(3.14) 

3m2b4  -  (15m2  -  3)b5  -  3b8  -  15b„  =  AA0m5  -  6AA2m2  -  4AA4m4 


(3.15) 

(m2  -  36) b6  +  bT  =  AA,  (3.16) 

21m2b6  -  (6m2  +  21)b7  -r  6b8  =  3AA2  (3.17) 

10m2b7  +  (15m2  +  10) b8  +  15b9  =  6AA2m2  -  SAAjm'  (3.18) 

3m2b3  *  (28m2  +  3)b<,  =  SAAom1  -  2AAjm6  (3.15) 


Partially  differentiating  equations  (3.1)  through  (3.4),  substituting 


into  equation  (2.39)  and  equating  the  coefficients  of  like  terms  yields: 
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(120m4  +  144m2  +  120)C„  +  (-144m2  -  240)C4  +  (-240m2  -  144)C2 

(3.20) 

+  120C3  +  72C4  +  120C5  =  8Aj  -  24A§m2  -  32AQAJ 
-4A2  -  12A(A2  +  4Ab0m2  -  2Ab2 

(120m1  +  480m2  +  840) C,  -  (240m2  +  240) C2  +  (-480m2  -  1680) C3 

(3.21) 

+  (-240m2  -  480) C4  -  240C5  +  840C,S  +  240C7  +  120C8 
=  SOAqAj  -  48A0A1m2  +  16A2  +  48Agm2  +  16A2  +  48A,A2 

-  2Ab4  +  4Ab2  +  4Ab[m2  -  4Ab0m2 

(240m1  -  240m2)  C,  +  (840m1  +  480m2  +  120)C2  -  240m2C3  (3.22) 

+  (-480m2  -  240) C4  +  (-1680m2  -  480) C5  +  120C7 
+  240C8  -  840C9  =  240AoA|m2  16A0Ajm4  -  32AoA2 
-  256A0A,m2  -  8OA20m2  +  144A|m4  -  112A,A2 
+  192A4A2m2  -  48A2  -  32A2m2  -  4Ab5  +  8Ab2m2  -  4Ab0mJ 
(120m1  +  1008m2  -  3024)C3  -  (240m2  +  504)  C4  +  120C5  (3.23) 

+  (-1008m2  -  6048)  C6  +  (-240m2  -  1008)  C7  -  240CS 
+  3O24C,0  +  504Cn  +  120C,2  =  -24Ajj 

-  24A20m2  +  96AoA,m2  -  48A2  -  24A)m2  -  72A,A2 

-  2Ab7  -r  4Ab,  -  2Ab2  4  4A  ,3m2  -  4Ab5m2 

(800m1  +  1680)  C3  +  (840m1  +  1600m2  +  840)C4  +  ,','380m2  +  800)C5 

(3.24) 

-  1680m2C6  +  (-1600m2  -  1680) C-  +  (-1680m2  -  1600) C„ 
-1680C,  +  840C]  [  +  800C12  +  840C1?  =  -640AaA1m2 

+  256A0A,m1  +  64AgA2  +  512A0A2m2  +  400AtA2  -  832A1A,m2 
+  48Agm2  -  144 A gm1  -  16A2m2  +  16A2m4  -  4Abt!  +  8Ab5 
+  8Ab4m2  -  8Ab2m2  -  4Ab4m4 
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120m<CJ+  (504m‘  ♦  240m*)Cj  ♦  ,3024m<  +  1008m1  *  120)CS  ,3.25) 

-  240m  C7  +  (-1008m2  -  240)C8  +  (-6048m2  -  1008)C3 

-  120C12  +  504C,3  +  3024Cn  -  72A^m4  -  120Agm6 

-  384AaA1m4  -  OBA^m6  -  192A0A2m2  +  76SApA2m4 

+  4S0A^  -552A jA2m4  -  72A2m4  +  24A|  -  24A22m2 

-  6Abq  12Ab5m2  -  6Ab2mJ 


(120m4  +  1728m2  +  7920)  Ce  + 

+  (-1728m?  -  15,  840) C 


(240m2  t  864)  C-  +  120C„ 
4  8 

!9  +  (-240m2  -  1728)Cn 


(3.26) 


-  240C12  =  -  48AgAj  -  48AoAim2  +  96A2  +  48A2m2 

-  48A,A:  -  4Ab-  -  2Ab4  +  4Ab6m2  -  4Ab3m3 

(1680m4  +  6048m2)  C6  +  (840m4  +  3360m2  +  3024)Cr  (3.27) 

+  (1680m2  +  *680)  C„  +  840C3  -  6048m2CIO 
+  (-336°m2  -  6048) Cn  +  (-1680m2  -  3360)C12 

-  1680C,3  =  336A0Aim2  -  272A0A,m1  -  96A„/\2 

-  256A0A2m‘  -  144A2m2  +  112A?m'  -  336A,A2 

-  1088A<A2m2  -  144A22  -  8Ab7m2  +  8Ab8  -  4Ab5 

-  8Ab4m2  _  4Ab3m4 

840m4Cg  +  ,1680m4  +  1680m2) Cr  +  (3024m4  ,  3360m2  +  840)  C8 

+  ,6048m!  *  168<»C.  -  *  (-3360m*  -  1680, 

*  (-6048m*  -  3360) C,;[  -  6048C,,  =  SSOA^.m1 

-  144A|jA,m6  +  64AoAzm'  -  IOSAjA;,,,1  -  1280A,A,m2 

+  1872A'A*mi  -  *6A*m*  -  96A*m«  +  48a|  +  48A(m* 

12Abg  +  12Abgm2  -  12Ab5m2  -  6Ab4m' 


t 
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120m4C7  +  (864m4  +  240m2)C8  +  (7920m4  +  1728m2  +  120)C9  (3.29) 

-  240m2C12  +  (-l/2Sm2  -  240) C13  +  (-15,840m2  -  1728)  C,4 
=  176A3A,mS  +  SOAoAjm8  +  IGOA^m4  -  512AoA2m6 

-  624A,A2m4  +  512A1A2m6  +  64A2m6  -  112A|m2  +  352A2m4 
+  IBAbgm2  -  dAb5m4 

(120m1  +  2640m2  +  17,16O)C10  +  (240m2  +  1320)0,!  (3.30) 

+  120C12  =  -SOA-  -  24A2m2  -  12A,A,  -  2Ab7 

-  4Ab6m2 

(2880m4  +  15,  840m2)  C10  +  (840m4  +  5760m?  +  7920)C„  (3.31) 

+  (1680m2  +  2880)  C12  +  840C13  =  96A2m2  -  128A2m4 
-f  48A,A2  -  448A,A2m2  +  48A2  -  4Ab8  -  8Ab7m2 

-  4Absm4 

3024m4C)0  +  (3528m4  +  6048m2)  C„  +  (3024m1  +  7056m2  +  3024)  C12 

(3.32) 

+  (6048m2  +  3528)  C13  +  3024CU  =  192A2m4  -  24A?m6 
+  672A,A2m2  -  1320A,A2m4  -  72A2  -  24A*m2  -  6Ab9 
••  l?Abgm2  -  6Ab7mJ 

840m4C,,  +  (2880m4  +  1680m2)  C!2  +  (7020m4  +  5760m2  +  840)  C,3 

(3.33) 

+  (15,840m2  2880)C14  =  16A]m6  +  80A?m8 

+  880A,A2m4  -  1024A,A2mG  +  16A2m2  -  352A2m4 

-  16Ab,,m2  -  8Absm4 

120m4C,2  +  (1320m4  +  240m2)C,3  +  (17,160m4  +  2640m2  +  120)C14 

(3.34) 

=  256A1A2mG  -  140A,A2m8  +  SSAjqn4  -  280A2me 

-  2QA2m8  -  10Ab,,m4 
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really  differentiating  equations  (3. 1)  through  (3. 5) ,  substituting  into 
equation  (2.42)  and  equating  '.ho  coefficients  of  like  terms  yields: 


<-2m?  -  2)d0  +  2d,  +  2dj  ^  A0b0 


-  12)d>  -  +  12d»  +  2*<  ■  -*AC0  +  3A0bj  -  5AA  t  AA 


(3.35} 


(~2m!) dj  +  (-12m?  -  2)d;  -  2d,  +  12ds  =  2ACom!  .  7A^, 


(3.36) 


(3.37) 


+  A0b2  +  3A2bo 

(-2^-  -  30)d3  -  2d,  +  30dt  +  2d,  «  -2AC,  +  4AC,  +  oA,bs  (2.38) 
-  HA0b,  +  7A0b0  +  3Ajb,  -  5A,b0 

M2mVi  +  (-12:r'!  '  -  12^  *  M4,  *  12d,  .  -6AC,  (3.39) 

+  6ACim2  +  3Ao*>4  “  5A0b2  -  21A0b,m2  +  l8Aobom2 
-  3A,b2  -  3A,b0m2  +  9A2b,  -  I5A2b0  +  8A0b, 

( -2m2) d4  +  (-30m2  -  2) d5  +  2d8  +  30d9  =  2AC2m2  -  4AC0m<  (3.40) 

+  A0b5  -  7A0b2m2  +  llA0b0m4  +  3A2b*  -  13A2b„m2 
(-2m2  -  56) d6  -  2d7  +  56dI0  +  2dn  =  4AC,  -  2AC3  -  2AC0  (3.41) 

7A0b6  -  17A0b3  +  13A0b,  -  3A0b0  +  5A,b3  -  HA,b, 

+  7A,b0 

(-30m2) d8  +  (-12m2  -30)d?  -  12ds  +  30dn  +  12d.  -  12AC2  -  6AC4 

-  8AC,m2  -  2AC0m2  +  10AC3m2  +  5A0br  -  3A0b4 

-  35A0bsm2  -  9  A0b2  -  llA0b0m2  +  46A0b,m2  -  A,b4 

-  35A,b2  -  17A,b,m2  +  10A,b0m2  +  l5A2b3  -  33A2b,  4  2lA2b0 


(3.42) 
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( '12m2) df  +  (-30m2  -  12)  dg  -  30d9  +  12d12  +  30d13  =  8AC2m2  (3.43) 

-  10AC5  +  6AC4m2  -  12ACim4  +  2AC0m4  +  3A0b8 

+  HA0b5  -  21A0b)m2  +  33A0bjm4  -  13A0b0m4  +  2A0b2m2 

-  7Ajb5  +  9A.(b2m2  +  3A,b0m4  +  9A2b4  -  7A2b2  -  39A2btm2 
+  38A2b0m2 

(-2m2)d8  +  (-56m2  -  2)dg  +  2d13  +  56d14  -  2AC5m2  -  4AC2m4  (3.44) 
+  2AC0m6  +  A0b9  +  llA0b2m4  -  5A0b0m6  -  7A0b5m2 
+  3A2b5  -  13A2b2ir.2  +  17A2b0m4 

(-2m2  -  90) d10  -  2dn  +  90d15  +  2die  =  4AC3  -  2AC6  -  2AC4  (3.45) 

-  23A0bg  +  19A0b3  -  5A0bj  +  7A4bg  -  17A4b3  +  13Ajb4 

-  3Ajb0 

(-56m2)d10  +  (-12m2  -  56) dlt  -  12d12  +  56die  +  12d17  =  12AC4 

(3.46) 

-  6AC7  -  16AC3m2  +  2AC1m2  -  6AC2  +  14AC6m2 

-  25A0b1m2  -  9A0b7  -  49A0b6m2  -,3A0b4  -»  5A0b2 
+  74A0b3m2  +  A4b7  +  gA^  -  31Ajb3m2  -  21A,b2 

-  7A]b0m2  +  SSAjbjin2  +  23A2b6  -  51A2b3  +  39A2b;  -  9A2b0 
(-39m2) dn  +  (-30m2  -  30)d12  -  30d13  +  30d17  +  30d18  =  20AC5 

(3.47) 

-  10AC8  -  10AC2m2  +  lOACjm4  +  10AC7m2 

-  20AC2m  +  5A0bg  -  35A0b7m2  -  25Agb3  +  5A0b2m° 
55A0b3m4  -  SSAob.m1*  +  3OA0b,,m2  -  5Ajb8  +  35A,b5 

-  5 A  1b4ni''  +  25A,b,m!  -  SA^m4  -  30A]b2m2  +  15A2b7 

-  25A2b4  -  65A2b3m2  +  5A2b2  -  25A2b0m2  +  90A2b4m2 


> 
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(-12m2)d12  +  (-56m2  -  12)d13  -  56d14  +  12d18  +  56d13  =  16AC5m2 
MAC 9  -  2AC2mJ  +  6ACsm?  -  12AC,mJ  +  6AC,mB  ^ 

+  19Aflb»  "  21AobSm2  +  33A0b4mJ  -  5A0b2m4  -  15A0b1m6 
MA0b5m2  -  llA,bg  +  21A jb5m2  -  9A,b2m4  -  A,b0me 

+  9A2b8  +  A2b5  -  39A2b.m2  +  51A2b,m4  -  2SA2b0m4 

+  22A2b2m2 

(-2m2) di3  +  (-90m2  -  2)d14  +  2dl3  +  9Od20  =  2AC9m2  -  4AC5m' 

2AC2m6  -  7A0b9m2  +  llA0b5m4  -  5A0b2m6  +  3A,b3  (  ’ 

-  13A2b5m2  +  17A2b2m;  -  7A2b0me 

(-2m  -  132 ) dj g  -  2d16  -  132d2,  +  2d22  =  4ACe  -  2AC10  (3.50) 

-  2AC.J  +  25A0be  -  7A0b:j  -  23A,b6  +  19A,b:?  -  5A,b, 
(-90m2)d15  -  (-12m2  -  90)d16  -  12dIT  +  90d22  -  12d,;J  -  l2AC: 

-  6ACn  -  24AC6m2  -  6AC ,  +  6ACjm2  -  18ACinm2  ('3’51) 

~  3A0b-  -  3A„b|  -  39A0b:im2  +  102Apbem2  +  3A,b7 

-  45A]b6m“  -  15A,b,  -  9A,b2  -  21A1b1m2  -  66A,b;im2 

-  69A,bG  +  57A2b:j  -  15A,b, 

(-56m2)  djg  -  ( -30nr  -  56)di:  -  30d18  .  56d,;j  .  30d21  (3.52) 

- 1 OAC  12;  +  20AC8  -  8AC;m2  -  10AC5  -  6AC,m2 
ISACjm1  ^  14ACnm2  -  28ACfim'  -  19A0b8 
13A„b5  -  '!Anb,nr  *  77Anb,-im1  -  57A0b,ml 
+  58AoMi2  <  29A,b8  -  19A,b7m2  -  49A,b3  •  21A,b2m2 
*  47Alb3m *  -  27A|blm1  -  2A]b,m2  -  43A2b7 
-  91A2b6m2  23A,bj  -  A2b2  -  51A,b,m2  ■  3 12A?().?m? 
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(-30m2)d47  +  (-56m2  -  30)d48  -  58d,3  +  30d24  +  56d25  (3.53) 

=  -14AC13  +  2 SAC 9  +  8AC8m2  -  18AC5m2 
+  6AC4m‘  +  10AC12m2  -  20AC7m4  +  10AC3m6 

-  41A0b9  +  21A0b5m2'  +  SSAofym1  -  27A0b4m4 

-  25A0b3m6  -  l4A0b8m2  +  55A4b9  -  TAjbgm2 

+  13A,bim4  *  15Aib2m‘  -  llA^m6  -  70A1b5m2 

-  17A?b8  -  e5A:brm2  -  llA2b5  -  9A2b2m2 
+  85A2b;im!  -  57A2b,mJ  +  74A2b4m2 

(-12m2) d18  +  (-90m2  -  12)d19  -  90d2n  -  12d26  +  90d2S  -  -18AC14 

(3.54) 

+  24AC:)m2  -  6AC5m’  +  6AC13m2  -  12AC8m4 

-  6ACjinG  -  33Aut8m4  i-  3A„b5m4  -  15A0b4mG 

-  30A9b9m2  +  SSAjbyin^  -  2SA]b:/n1  -  3A,b2ms 

1  I 

~  9A2b9  -  39A2b8mw  -  51A2b|in  -  15A2b2m 

-  21A"h1m6  +  6A2b5m2 

(-2m2)d19  +  (-132m2  -  2)d2n  +  2d26  +  132d27  =  ?ACum2  (3.55) 


-  4ACflm'  - 

2AC5m8  + 

llAob^m1  -  5A0b5mG 

-  13A2b9m2  ■+ 

•  ITAjbspi1 

-  7A2b2mh 

(-2m2  -  182) d2,  -  2d22 

•  lUSdjg  1 

' ,  2dj.j  =  4AC|0  -  2A'  8  (3. 

-  9Anb(i  •  25A,b 

s  -  7A,b, 

(-132m2!  d2l  *  (-12m2  - 

I'Alit  1-.o  •* 

1  !32dj.,  •  1  4(0  (3. 

12AC - 

32AC„im 

ii.U'  •  •  l'tA.C ,>!>'• 1  '  A()b; 

-  58Ajb{inr' 

- 

~ 1  , ,  h  (Mi*1  ■  9 1.\  |b,:m 

•  75A -I),;  - 

21 A  .3. 
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(••90m2)  d22  4  (-30m2  -  90) d23  -  30d24  +  90^  +  30d31  (3.58) 

=  20AC12  -  16AC14m2  -  10ACg  -  2AC7m5 
+  26ACdm4  -  36AC10m4  +  HA^bg  -  23A0b7m2 

-  79A0bgm4  -  43A]b8  +  21A1b5  +  7Atb4m2 

+  69Ajb5m4  -  49A]b3m4  +  26A5b7m2  +  41A2b; 

-  7A2b4  -  77A2b3m2  *  194A2bem2 

(-56m2) d23  +  (-56m2  -  56)d24  -  56d23  +  56d3i  +  56d32  (3.59) 

=  28ACjj  -  14 AC 3  -  14AC8m2  +  14ACrm4 

-  28ACjjm4  +  14ACgm6  +  21A0b9  +  7A0b8m2 

-  49A0b7m4  -  35A0bfimG  -  77A1b9  +  49Ajbsm2 

+  35Ajbjm4  -  7A1b1mJ  -  21A1b3m6  -  42A1b8m2 
+  7A2b8  +  7A2b5  -  35A2b5m2  119A2bsm4 

-  91A2b3m4  +  126A2b7m2 

<  -30m2)  d24  +  (-90m2  -  30)d,5  -  90d26  +  30d32  +  90d33  =  36ACU  (3.60) 
+  16AC)3m2  -  2SAC9m2  +  2AC8m4  -  20ACj2m4 
+  10AC7m6  +  37A0bsm2  -  19A0b8m4  -  25Ayb7m6  +  Ajb8m4 
+  35Ajb5m4  -  7Ajb4m6  -  110A1b3m2  -  £7A2b9  +  7A2b5m2 
+  85Aib7m4  -  49A2b4m1  -  35A2b3m6  +  58A2bsm2 
(-12m2)d25  +  (-132m2  -  12)d26  -  132n27  +  12d;i:,  +  132d34  (3.61) 

-  32ACi4m2  -  10AC8m4  -  12ACl3m4  +  6ACgin6 
+  llA0b9m4  -  15A0b8m6  -  33A5b9m4  t  7A1b5m6 
+  51A2bgm4  -  7A2b5m4  -  21A2b4m6  -  1°A2b9m2 


(-2m2)d26  +  (-182m2  -  2)d27  +  2d34  +  182d35  =  l?A2bs)m4  (3.62) 
-  SAobjm6  -  7A2b5m6  +  2AC9m6  -  4AG14m4 
(—2m  —  240)  d»j  —  2d29  —  — 9A4hg  —  2AC48  (3,  •  >3 ) 

(-182m5) d28  +  (-12m2  -  182) d29  -  12d30  =  14AC10m2  -  6AC44  (3.64/ 
+  5A4b7  -  49A1bgmi  -  27A2'os 

(-132m?)d29  +  (-30m2  -  132)d30  -3Gd31  =  2ACnm2  -  10AC12  (3.65) 

+  34ACI0mI  -  19Ajb8  -  7Ab?m2  -  71A,b6m4 


-  13A2b7  ~  103A2b6m2 

(-90m2)  dl0  +  (-56m2  -  S0)d31  -  56d32  =  -14AC,3  -  10AC12m2  (3.66) 

-  22ACnm4  +  18ACi0m6  +  33Aib9  +  354,b8m2 
+  A2b8  -  29A]b7m4  -  3-lAjbgni6  -  61A2b7m2 

-  125A2bem4 


(-56m2)d31  i-  (-90m2  -  56)d32  -  90d33  =  -18AC14  -  22AC13m2  (3.67) 

+  10AC1?m4  +  MACijm*  +  77A,b9m2  +  13A,b8m4 

-  r7Ajb7m6  »-  15A2b9  -  19A2b8m2  -  83A2b7m4 

-  49A2bsm6 

( -30m2) d32  +  (-132m?  -  30)d33  -  132d34  =  -34AC14m2  (3.68) 

-  2AC13m‘  +  10AC12m6  +  55A]b9m4 

-  3A,b8m6  +  23A2b9m2  -  41A2bsm4  -  35A2b7m6 

(-12m2) d33  +  (-182m2  -  12)d34  -  182  d35  =  -14AC14m4  (3.69) 

+  6AC13mG  +  llAjb9m6  +  A2bKm4  -  21A2b8m6 
(-2m2)  d34  +  (-240m2  -  2)d35  =  2AC14m6  -  7A2b9m6 


(3.70) 


CHAPTER  IV 


COMPUTER  ANALYSIS  AND  NUMERICAL  EXAMPLES 


Equations  (3.7)  through  (3.9)  for  the  first-order  approximation  of  y 
are  written  in  matrix  form  and  designated  as  the  W-matrix  with  the  following 
nonzero  elements: 


Wj  i  =  5m4  +  2m2  +  1 

Wj  2  =  -2m2  -  2 

Wt  j  =  -10m2  -  2 

w  A2m2 

W14-  12 

W2  j  =  10m4  +  6m2 

W2  3  =  105m4  +  20m2  +  3 

wj<  -  -  AY 

W3  2  ~  5m4  f  12m2  +  15 

W3  3  =  10m2  +  6 

w  -  A-2P*-- 

*34  12 

wM  w1>2  W1>3 

Ao 

I - 

W2(1  W2)2  W2>3 

A! 

1! 

W,  4 
'  » 

'  3 

CO 

£ 

_a2_ 

_W3,4_ 

W-Matrix 


Equations  (3.10)  through  (3.19)  for  the  first-order  approximation  of 
W  are  written  in  matrix  form  and  designated  as  the  X-matrix  with  the  following 


nonzero  elements: 
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X,, 

=  m2  +  3 

X 

II 

-3 

X,  3 

=  -1 

X,  „  = 

=  AA0 

X*i  2 

=  m2  +  10 

X23  = 

1 

X|4 

=  -10 

X25  = 

-1 

X2  11 

=  AAj  -  2AA0 

CM 

CO 

X 

=  3m2 

X;l  3  = 

6m2  +  3 

X3  5 

=  -3 

Xs  6  = 

-6 

X3  11 

=  3AA2  -  2AA0m2 

-  2AA,m2 

X44 

=  m2  +  21 

X45  = 

1 

X4  7 

=  -21 

x48  - 

-1 

X4  11 

=  AAG  ~  2AAj 

X54 

-  2 
=  Dtn 

X55  = 

3m2  +  5 

x5G 

=  3 

X58  = 

-5 

X5  9 

=  -3 

X5,1  = 

AA0m2  + 

XG  5 

=  3m2 

Xg  g  = 

15m2  +  3 

xG9 

-  -3 

Xg  io  ~ 

-15 

XG  11 

=  AA0m!  +  4AA1m4 

-  6AA.2m2 

X7  7  = 

=  m2  +  36 

X7  3  = 

1 

X7  11 

-  AA] 

Xe  7  = 

=  21m2 

V  __ 

-“8  8  ~ 

6m2  +  21 

X89  =  6 


X8  11  -  3AA  ; 
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X-Matrix 


Equations  (3.20)  through  (3.34)  for  the  second-order  approximation  of 
i l>  are  written  in  matrix  form  and  designated  as  the  Y-matrix  with  the  following 
nonzero  elements: 

Y,  i  =  120m4  +  144m2  '  120  Y,  2  =  -144m2  -  240 

Y,  3  =  -240m2  -  144  Yj  4  =  120 

Yj  5  =  72  Y,  6  =  120 


Yj  16  =  8Ag  -  24A§m2  -  32A,>A4 

-  4A?  -  12A4A2  +  4Ab0m2  -  i 

Y2  2  =  120m4  +  480m2  +  840 

Y2  3  -  240m2  +  240 

Y2  4  =  -480m2  -  1680 

Y  =  -240m2  -  480 

2  5 

Y2  g  =  -240 

Y2  7  =  840 

Y, 8  =  240 

Y2  9  =  120 

Y2  ig  =  80AqAi  -  48AoA,m2  +  16A2  +  48A2m2  +  16A^  +  48A,A 

-  2Ab4  +  4Ab2  +  4Ab!m2  ■ 

-  4Ab0m2 

Y3  2  =  240m4  +  240m2 

Y3  3  =  840m4  +  480m2  +  120 

Y3  4  -  -240m2 

Y3  5  =  -480m2  -  240 

Y3  g  =  -1680m2  -  480 

Y3  8  =  120 

Y,  ,  =  240 

Y3  jo  -  840 

Ys  16  =  240A0A1m2  +  leAgAjm4 

+■  32A()A2  -  256AoA2m2 

-  80Agm2  +  l44A5m!  -  112A4A2  +  192A4A2m2  -  48A2 

J'  32A4m2  -  4Ab3  +  8Ab2m 

2  -  4Ab0m4 

Y4  4  =  120m4  +  100 8m2  -r  3024 

Y4  5  =  240m2  +  504 

Y,  4  =  120 

Y4  7  =  -1008m2  -  6048 

Y,  j  =  -240m2  -1008 

Yj  3  =  -240 

Y,  ,  -  3024 

Y,  12  =  504 

Y  |  =  120 

Yj  jg  -  -24Ag  -  24Aom2  +  96A0A,m2  -  48A2  -  24A2m2  -  72A4A2 


-  2Ab7  +  4Ab|  -  2Ab2  +  4Ab3m2  -  4Ab4m2 
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Y5  4  =  800m4  +  1680 
Ys  6  =  1680m2  +  800 
Y5  8  =  -1600m2  -  1680 
Y5  10  =  -1680 
Y5  13  =  800 


Ys  5  =  840m4  +  1600m2  +  840 

Y5T  =  -1680m2 

Ys  9  =  -1680m2  -  1600 

Ys  12  -  840 

Y5  14  =  840 


v5  „  =  256A^,m*  -  64OA0A,m2  4  64A.A,  4  SKA^m*  4  400A,A,X 
-  S32A,A;>m2  4  48ASm2  -  144A2m4  -  164fm! 

4  leAjm*  -  4Ab,  4  SAbs  4  SAb.m*  -  8Ab,m2  -  4Ab,m‘ 


Ye  i  =  120m4 

Ye  6  =  3024m4  +  1008m2  +  120 
Yg  g  =  -1008m2  -  240 
Y6  13  =  120 
Yg  15  =  3024 


Yg  5  =  504m4  -f  240m2 
Ye  8  =  -240m2 
y6  jo  =  -6048m2  -  1008 
Ye  I4  =  504 


Yg  ,g  72A20m4  -  120A^6  -  384A0A,m4  -  o6A0A1mB  -  192A0A2m2 
+  768A0A2m4  +  480A1A2m2  -  552A,A2m4  -  72A2m4  +  24A| 
24A2m2  -  6Ab9  +  12Ab5m2  -  6Ab2m‘ 


Y7  7  -  120m4  +  1728m2  +  7920  Y7  8  =  240m2  +  864 

Y? 9  ]2°  Y7  n  -  -1728m2  -  15,840 

Y7  »  =  -240m2  -  1728  Y?  „  =  _240 

Y7  ,s  -  -48A0A,  -  4 8AqA jin2  +  96A?  ;  48A?m"  >  48A,A2  +  4Ab 

-  2Ab4  +  4Abem2  -  4Abgm2 


Y8  7  -  16?  Jm4  +  6048m2  Y8  „  =  840m4  +  3360m2  +  3024 

Y8  9  =  1680m2  +  1680  Y8  i0  =  840 

Y?  ,4  =  -6048m2  Y8  12  =  -3360m2  -  6048 

Yg  ,3  =  -1680m2  -3360  Y8  14  =  -1680 

*8  ;s  =  3f  PAgAjm2  -  272AoAjm4  -  96AgA2  -  256A0A2m2  -  144A^m2 
+  lo.2A2m4  -  336A4A2  +  1088A4A2m2  -  144A2  +  8Ab7m2 
+  8Ab8  -  4Ab5  -  8Ab4m2  -  4Ab3mJ 

Y9  7  =  840m4  Y9  8  =  1680m4  +  1680m2 

Y9  g  =  3024m4  +  3360m2 +  840  Y9  10  =  6048m2  +  1630 

Y9  12  =  -1680m2  Y9  13  =  -3360m2  -  1680 

Yg  14  =  -6048m2  -  3360  Y9  15  =  -6048 

Yg  ls  =  560AgA,m4  -  144A0A1m6  +  64A()A2m2  -  768AgA2m4 

-  1280A4A2m2  +  1872A4A2m4  -  16A^mJ  -  96A^m6 
+  48A2  +  48A2m2  +  12Ab9  +  12Ab8m2  -  12Ab5m2 

-  6Ab4m4 

Yjo  8  =  120m1  Y10  9  =  864m 1  +  240m2 

Y10  io  =  7920m1  +  1728m2+  120  Yi0  i3  =  -240m2 

Y10  14  =  -1728m2  -  240  Y10  15  =  -15,840m2  -  1728 

Yjo  i6  =  176A0Aim6  +  80AgA,m3  +  leOAgAjm1  -  512AgA2mG 

-  624A,A2m4  +  512A!A2m6  +  64A2m6  -  112Aj;m2 
+  352A2m4  +  16Ab9m2  -  SAbjm1 
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Yn  n  =  120m4  +  2640m2  +  17, 160  Y„  12  =  240m2  +  1320 
Y«  13  -  120 

Yu  16  =  -60A2  -  24A2m2  -  12AjA2  -  2Ab7  -  4Abem2 

Y12  jj  =  2880m4  +  15,840m2  Y12  ,2  =  840m4  +  5760m2  +  7920 

Y12  13  =  1680m2  +  2880  Y12  14  =  840 

Y]2  16=  96A2m2  -  128A2m4  +  48A,A2  -  448Aji42m?  +  48A?  -  4Ab8 


-  8Ab;  -  4Abem4 

Y,3 

ii  =  3024m4 

Y 13  )2  =  3528m4  +  6048m 

13  =  3024m4  +  7056m2  +3024 

Y,3  ,i  =  6048m2  +  3528 

Y.3 

15  =  3024 

Y,3  i6=  192A2m4  -  24A2m6  +  672AjA2m2  -  1320AiA2m4  -  72A| 

-  24A2m2  -  6Ab9  -  j.2Ab8m2  -  6Ab7m4 

Yu  12  =  840m4  Yu  1:J  =  2880m4  +  1680m2 

Y,4  ,4  =  7920m4  +  5780m2-*-  840  Yi4  15  =  15,840m2  +  2880 

Yj4  )6=  16A2m6  +  80A2m8  +  880AjA2m4  -  1024A.A2ms  *  16Aj;m2 

-  352A2m4  -  16Abgm2  -  8Ab8m' 

Yi5  ,3=  ^Om4  Y 1 5  ,4  -  1320m’  +  240m2 

Y  i5  15  =  17, 160m4  +  2640m2  +120 

Yi5  ,6  =  256AjA2m8  -  140A,A2mB  +  88A2m’  -  280A22mb  -  20Ami8 

-  10Ab9m4 


Y-Matrix 


Equations  (3.35)  through  (3.70)  for  the  second-order  approximation  of 
,v  are  written  in  matrix  form  and  designated  as  the  Z-matrix  with  the  following 
nonzero  elements: 

Z|  |  -  -2m2  -  2  Zj  2  =  2 


Z,  ,  - 


Zi  37  =  A0b{) 


2 
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Z2  2  -  -2m’  -  12 
Z2  4  =  12 

Z2  37  =  -2AC0  +  3A0bj 


Z2  3  -  -2 
Z2  5  =  2 
5A0b0  -f  Ajb0 


z3  2  -  -2m2 

Z3  3  =  -12m2  -  2 

Z3  5  =  2 

Z3  e  =  12 

Z3  37  -  2AC0m2  -  7A0b0m2  +  A0b2  +  3A2b0 

Z4  4  =  -2m2  -30 

Z45 - 2 

Z4  7  =  30 

Z4  8  =  2 

2-4  37  =  4AC0  -  2AC]  4-  5A0b3 

-  llA0b,  +  7A0b0  +  3A,b] 

Z5  4  =  -12m2 

Z5  5  =  -12m2  -  12 

zs  e  =  -12 

Z5  £  =  12 

Z5  9  -  12 

Z5  37  =  6AC,m2  +  3A0b4  -  6AC 

2  -  5A0b2  -  21A0b1m2  +  : 

-  3A,b2  -  SA^m2  +  9A2b,  -  15A2b0  +  8A0b2 

?6  5  =  -2m? 

ZG  6  =  -30m2  -  2 

^69  =  2 

Z6  i0  =  30 

"6  37  -  2AC2m2  -  4AC„m4  +  A0b5  -  TA^m2  +  llA0b0m4 

-  13A2b0m2 

:7  7  =  -2m2  -  56 

Z7  8  =  -2 

'7  11  “  56 

Z7  t2  -  2 

7  37  =  4AC,  -  2AC3  -  2AC0  + 

7A0Og  -  17Agbj  +  13A„b, 

+  5A]B3  -  HA,b,  ►  7A,b0 

kiDo 


k2o2 


-  3Anb 


ouo 


% 
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Z8  7  =  -30m2 

Zg  9  =  -12 
ZR  13  -  12 


Z8  g  =  -12m1  -  30 
Z3  12  =  30 


Zg  37  =  12AC2  -  6AC4  -  8ACjm2  -  2AC3m2  ■+  10AC3rn2  +  5A0b7 

-  3A0b4  -  35Arb3m2  -  9Af,b2  -  llA0b0m2  +  46A0b1m2 

-  Ajb4  +  15Ajbj  -  ITAjbjm2  +  lOAjbgm2  +  15A.2b3 

-  33A2bt  +  21A2b0 


Z9  8  “ 

-12m2 

Z9  s  =  -3 'Jin2  « 

» 22 

Z9  10  = 

-30 

7  —  19 

13  “ 

Z9  14  = 

30 

Z#jT  =  SAC2ir.2  -  lOACj  +  6AC4m2  -  12ACjm4  +  2AC  m4  +  3A0b8 
llA„b,  -  21A,-,b1m2  +•  33A0bjm4  -  I3A0b0m4  +  2A0b2m2 

-  7  Ajb5  +  9Ajb2m2  +  3Ajb0m4  +  9A2b4  -  7A2b2  -  39A2b1m2 
+  38A2b3m2 

Zjo  9  =  _2m  Z,0  u  =  -56m‘  -  2 

ZlO  14  ~  2  Z13  15  =  r,S 

Zjo  37  =  2ACsm2  -  4AC\m!  +  2AC0m6  +  A0b9  +  31A0b2m4  -  5A0b0mG 

-  7A0b5m2  +  3A2b5  -  13A2b2m2  4  17A2bcm4 

11  ~  — 2mw  —  90  Zn  jo  —  -2 

Z'ii  16  ~  90  Zjj  j(  =  2 

Zn  37  =  4AC3  -  2ACe  -  2AC?  -  23A0bs  +  19A0b5  -  5A0bj  +  7Ajbe 

-  17A ;b3  +  13A,b]  -  3Ajb0 
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Z12  ii  -  -56m2  Zj2  i2  =  -12m2  -  56 

Z12  13  =  -12  Zj2  1T  =  56 

Zi?  ig  =  12 

Zl2  37  =  12AC4  -  6AC7  -  16AC3m2  +  2ACjm2  -  6AC2  +  MACgm2 

-  25A0b1mi  -  9A0b7  -  49A0b6m2  -  3A0b4  +  5A0b2 
+  74A0b3m2  +  A,b7  +  9Ajb4  -  SlAjbgin2  -  21A,b2 

-  7A1b0m2  +  SSAjbjm2  +  21A2be  -  51A2b3  +  39AJ.3!  -  9A2b0 

2 

Zj3  12  =  -30m  Z53  13  =  -30m2  -  30 

Zi3  14  =  "30  Z13  18  =  30 

Z13  19  =  30 

Z 1 3  37  =  20AC5  -  lOACj  -  10AC2m2  +  lOACjm4  +  10AC7m2 

-  20ACsm4  +  5A0b8  -  35A0b7m2  -  25A0b5  +  5A0b2m2 

-  55A0n3m"  -  35A0b,m4  +  3OA0b4m2  -  5A,b8  +  35A,b-5 

-  5A,b4m2  +  25A]b1m4  -  5A]b0m4  -  SOA^m2  +  15A2b7 

-  25A2b4  -  65A2b3m2  +  5A2b2  -  25A2b0m2  +  90A2b,m2 


’ll  13  =  -12m2 

zi4  14  =  -56m2  -  12 

14  15  “  -30 

Z14  Jg  =  12 

14  20  ~  06 

zl4  37  ~  16AC 5m2  -  14ACg  -  2AC2m4  +  6AC8m2  -  12AC4m4 
+  6ACjm6  +  l?A0bs  -  21A0b8m2  +  3A0b4m4  -  5A0b2m4 
15A0b1m*5  -  14A0b5m2  -  HAjb9  +  21A]b5in2  -  OA^m4 

-  Ajbom6  4  9A2b8  +  A2b5  -  39A2b4m2  51A2b,m4 

-  23A2b0m4  +  22A2b2m2 
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Zi5  14  =  -2m2  Z15  15  =  -90m2  -  2 

zlS  20  =  2  Z15  21  •-  90 

zi5  37  =  2AC9m2  -  4ACsm‘t  +  2AC2m6  -  7A0b9m2  +  llA0b5m4 

-  5A0b2m6  +  3A2b9  -  13A2b5m2  +  17A2b2m4  -  7A2b0m6 

zl6  16  =  -2m2  -  132  ZI6  17  =  -2 

z16  22  =  132  Z!6  23  =  2 

zi6  37  =  4AC6  -  2AC,0  -  2AC3  +  25A0b6  -  7A0b3  -  23AtbG 
+  19Ajb3  -  5A]bj 

Z17  16  =  -90m2  Z17  17  =  -12m2  -  90 

Z17  18  —  -12  Z17  23  =  90 

Z17  24  ~  12 

Zn  37  =  12AC7  -  6ACn  -  24AC6m2  -  6AC4  +  6AC3m2  +  18AC10m2 
+  3A0b7  +  3A0b4  -  39A0b3m2  +  102A0bGm2  +  3A4b7  -  45A1b6m2 

-  15A,b4  +  9A4b2  -  21A,b1m2  +  OGA^m2  -  69A2b8  +  57A2b3 

-  15A2b4 

Zt8  I7  =  -56m2  Zlg  18  =  -30m2  -  56 

Z18  19  =  -30  zl8  24  =  56 

Z18  25  =  30 

Z18  37  =  20AC8  -  10AC,2  -  8AC7m2  -  10AC5  -  6AC4m2  +  18AC3m4 
+  14ACj]m2  -  28ACem4  -  19A0b8  +  13Anb5  -  9A0b4m2 
+  7/Agbgm  -  57A0b3iv.4  +  58A9b7m2  +  29A7b8  -  19Ajb7m2 
-  49Ajb5  +  21A,b2m2  +  47A,b3m4  -  27A1b1m4  -  2A1b4m2 
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-  43A2b7  -  91A2b6m2  +  23A2b4  -  A2b2  -  51A2b,m2  +  142A2b3n2 

=  "30m2  Zl3  19  =  -56m2  -30 

2  ia  ?o  ~  -56 

Z 19  26  -  56 


Z19  25  -  30 


37  =  28AC9  -  14AC13  +  8AC8m2  -  18AC5m2  +  6AC4m4 
+  10AC.2m2  -  20AC7m4  +  10AC3m6  -  41Aob9  +  21A0b5m2 

-  55A0b7mJ  -  27A0b4m4  -  25A0^m6  +  14A0b8m2  +  55A,b9 
7A,b8Ki‘  -  13A]b4m4  v  15A,b2m4  -  llA,b,mfi  -  70A,b5m2 

-  17A2b8  -  65A2b7m2  -  llA2b5  -  9A2b2m2  +  S5A2b3m4 

-  57A2b1m4  -  74A2b4m2 


Z2q  59  =  -12m2 
^20  21  ~  -90 
Z20  27  -  90 


z2o  20  _  -90m2  -  22 


Z20  26  -  12 


z20  37  -  24AC9m2  -  18AC) 4  -  6AC5m‘  +  6AC13m2  -  12ACsm4 

-  6AC4I116  +  33A0b8m4  +  3A0bsm4  -  15A0b,m6  -  3OA0bym2 

-  33Ajb9m2  +  28A1bim4  +  SA^m6  +  9A2b9  -  39A2b8m2 
51A2b4m  -  ISAjbjm1  -  21A2b1m6  +  6A2bsm2 


z2i  20  =  -2m2 


z2i  2i  -  -132m2  -  2 


Zjj  27  —  2 


Z21  28  “  132 


z2i  37  =  2ACum2  -  4ACsm4  1  2AC5m6  +  llA0b9m4  -  5A0b5mfi 
-  13A2b9m'-'  +  17A2b5m4  -  7A2b2m6 


z22  22  =  -2m2  -  182 


Z22  23  ~  “2 


I 


Z90  9Q  —  182 


Z22  30-2 


Z22  37  “  4AOio  “  2ACg  -  9A0bg  +  25A1bg  -  TAjbg 


Zio  m  —  -132m2 


ZM  -  -12 


Z23  23  =  -12m2  -  132 


Z2330  “  132 


Z23  31  -  12 


Z23  32  —  12ACjj  —  32ACjom2  —  6AC7  +  lOACgm  +  Agb7  —  53Agbgm 
-  9Ajb7  +  7Ajb4  -  SSAjbgm2  +  94A1bgm2  +  75A2b6  -  21A2b3 


Zo.  =  -90m2 


Z04  94  -  —30m.  —  90 


Z2J  25  -  -3° 


Z21  31  -  90 


Z?4  32  =  30 

Z21  37  =  20ACl2  -  16ACnm2  -  10AC8  -  2AC7m2  +  26ACGm4 
-  36AC10m4  -1-  llA0b8  -  23A0b7m2  -  79A0bem4  -  43A1b8 
+  21A1b5  +  7A1b1m2  +  egAjbgm  -  lOAjbgm4  -r  26A)b7m2 
+  41A2b7  -  7A2b[  -  77A2b:)m2  +  194A2b6m2 


Zo=  01  —  —  56m 


Z25  25  —  -5bm  —  56 


Zo?  ~  56 


^25  33  “  36 


Z25  :!7  =  28AC13  -  14AC(,  -  14ACgm2  +  14AC7m4  -  28ACnm4 
3  14ACeme  +  21A0b9  +  7Anb8m2  -  49A0b7m4  -  35A0bGm6 
-  77A1b<)  +  49A1b5m2  +  35Aib7m'  -  7A)b4m4  -  21A]b3m6 


-  42A!b8m2  +  7A2b8  +  7A2bs  -  SSA^jm"  +  119A2b6m4 

-  91A2b:jin4  +  126A2b7m2 


> 
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^26  25  30m"  Z2e  26  -  -90m2  -  30 

Z26  27  ^  -90  Z2g  33  =  30 

Z26  34  ~  9° 

z2e  37  =  36AC14  +  16AC13m2  -  26AC9m2  +  2ACsm4  -  20AC12m4 

-  10AC7mG  +  37A0b9m2  -  lSA,^4  -  25A0b7me  +  A^m4 

+  35A,b5m4  -  7A,b4me  -  110A,b9m2  -  27A2b9  +  7A2b5in2 
+  85A2b7m4  -  49A2b4m4  -  35A2b3m6  +58A2b8m2 

Z27  26  -  -12m2  Z27  27  =  -132m2  -  12 

Z2  7  28  =  -132  Z27  34  -  12 

Z27  35  =  132 

Z27  37  =  32AC14m2  -  10ACgm4  -  12AC13m4  +  6AC8m6  +  !.lA0b9m4 

-  15A0b8m6  -  33Ajb9m4  +  7A3b5m6  +  51A2b8m4  -  7A2b5m4 

-  21A2b4m6  -  10A2b9m2 

Z28  27  =  -2m2  Z28  28  =  -182m2  -  2 

Z28  35  =  2  Z28  3(3  =  182 

z28  37  =  17A2b9m4  -  5A0bt.m6  -  7A2b5m6  +  2AC9rn6  -  4AC14m4 

Z29  29  =  -2m2  -  240  Z29  30  =  -2 

z29  37  =  “9A]b6  -  2ACi0 

z30  29  =  -182m2  Z30  30  =  _12mz  _  182 

Z30  31  =  “12 

z30  37  =  14AC10m2  -  6ACjj  +  5A3b7  -  49A3b6m2  -  27A2b6 
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Zj,  30  —  -132m2  Z34  oj  —  -30m2  -  132 

ZS1  S2  “  -30 

ZS1  37  =  2ACnm2  -  10AC12  -  34AC10m4  +  WAjbg  -  TAjfym2 

-  TlAjbgm4  -  13A2b7  -  103A2b6m2 

Z32  37  =  —90m2  Z32  32  —  —56m2  —  90 

Z32  33  =  -56 

Z32  37  =  22AC77m  —  MAC  13  —  10AC72m  +  18AC7Qm®  +  33A]Dg 
+  SSAjbgm2  +  A2bg  -  29Ajb7m4  -  SlAjbgm6  -  61A2b7m2 

-  125A2b6m4 

z33  32  =  -56m2  Z33  33  =  -90m2  -  56 

z33  34  =  "SO 

Z33  37  =  10AC12m4  -  18AC14  -  22AC13m:  +  14ACnm6  +  TTAjbgm2 
-r  13Ajbgm4  -  MAjfym6  +  1  )A2b9  -  lSA^m"1  -  83A2b7m4 

-  49A2b6m6 

Z34  33  =  -30m  Z34  34  =  -132m2  -  q0 

Z34  35  =  “132 

Z34  37  =  10ACj2m  -  34ACj4m2  -  2ACj3m4  +  oSAjbgiii4  -  SAjbgm® 
+  23A2b9m2  -  41A2b8m4  -  35A2b7m6 

Z35  34  =  -12m2  Z35  35  -  -182m"  1. 

Z35  36  =  “182 

Z35  37  -  6ACl3me  -  l4AC14in4  +  llA^m6  4  A2bsm4  -  2L-3bPm6 


^36  35  ~  -2m2  Zjg  36  -  -240m2  -  2 

Zjg  37  =  2ACj4m6  -  7A2bam 

All  elements  in  the  W-,  X-,  Y-,  and  Z-matrices  which  are  not  other¬ 
wise  defined  are  zero. 

The  total  volume  rate  of  flow  through  the  pipe,  Q^,  is  given  by 
Qt  =  W0a2(Q0  +  KQj  +  where 


Qo 


4/ 


x=0 


y=— n/  1  -x2 
m 

/  W.j  dy  dx 

y=o 


x=l  m 


Q.  -  4  /  / 

x=0  y=0 


Wj  dy  dx 


Q2 


4 / 


x=0 


y-— n/  1  -x2 
m 

f  Wj  dy  dx 
y=0 


(4.1) 


(4.2) 


(4.3) 


Since  the  equations  are  in  nondimensional  form,  the  limits  in  the  x-direction 
are  taken  to  be  from  x  =  0  to  x  =  1  and  the  limits  in  the  y-direction  are  from 
v  =  0  to  the  boundary  of  the  cross  section. 

From  equation  (3.3)  it  is  noted  that  Wj  is  a  function  odd  in  x  and  even 
in  y  integrated  over  an  area  symmetric  about  the  origin.  Therefore,  Q!  =  0 
and  the  first  approximation  of  Whave  no  effect  on  the  flow  rate. 

Since  Qj  =  0,  the  equation  for  the  total  rate  of  flow  through  the  pipe  is 


reduced  '0: 
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Qt  =  W0a2[Q0  -  Q2K2]  . 

The  product  Q0W0a2  is  the  rate  of  flow  through,  a  straight  p;oe.  Therefore,  the 
total  flow  rate  is  expressed  as: 


Qt  ~  QoW0a2 


(4.4) 


The  bracketed  terms  represent  the  reduction  in  flow  rate  due  to  the  curvature 
of  the  pipe  and  are  designated  as: 


(4.5) 


where  F  /F  is  the  ^tio  of  the  flux  through  a  curved  pipe  to  the  flux  through  a 
c  s 

straight  pipe,  both  having  the  same  cross  section,  length  and  inlet  pressure. 

Equations  (3. 1)  and  (3.5)  are  even  in  both  x  and  y,  and  to  simplify 
integration  are  written  as  follows: 

W0  =  A  -  Ax2  -  AmV  (4.8) 

W2  =  d0  +  (dj  -  d0) x2  +  (d2  -  d0m2)y2  +  (ds  -  dj)x4  -r  (d4  -  d2  -  d,m2)x2y2 

(4.7) 

+  (d5  -  d2m2)y4  +  (d6  --  d3)x6  +  (d7  -  d4  -  d3m2)x4y2 
+  (d8  -  d5  -  d4m2)x2y4  +  (d9  -  d5nr)yG  +  (d10  -  d6)x8 
+  (d„  -  d7  -  d6m2)x6y2  +  (d!2  -  d8  -  d7m2)xV  +  (d13  -  d3  ~  d8m2)x2y6 
+  (d!4  -  d«,m2)y8  +  (d!5  -  d10)x10  +  (d,6  -  dn  -  dlcm2)x8y2 
+  (d17  -  d!2  -  dnm2)xy  +  (d18  -  d,3  -  d,2m2)x4y6 
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+  (<*i9  -  d14  -  d13m‘)xy  +  (d20  -  d14m2)y10  +.  (d21  -  d15)x12 
+  (<*22  -  di6  -  di5m2)x10y2  +  (d23  -  d17  -  diem2)x8y4 
+  (d24  -  d18  -  d17m2)x6yG  +  (d25  -  d19  -  d18m2)x4y8 
+  (d26  -  d20  -  dl9m2)x2y10  +  (d27  -  d20m2)y12  +  (d28  -  d21)x14 

+  (d29  -  d22  -  d2lm2)x12y2  +  (d30  -  d23  -  d22m2)x1V 

+  (d3i  ~  d24  “  d23m:)xV  +  (d32  -  d25  -  d24m2)x6y8 
+  (d33  -  d26  ~  d25m2)xJy10  +  (d34  -  d27  -  d26m2)x2y12  +  (d35  -  d27m2)y14 

+  ( -d28x16)  +  (-d29  -  d28m2)xI4y2  +  (-d30  -  d29m2)xI2y4 

+  (-d31  -  d30m2)x10y6  +  (-d32  -  d31m2)x8y8  +  (-d33  -  d32ir.2)x6y:o 
+  (-d34  -  d33m2)x4y12  +  (~d35  -  d34m’)x2yt4  +  (-d35m2)y16 

The  integrated  portion  of  each  term  of  equations  (4.6)  and  (4.7)  is  represented 

by: 

y-—  m'  1  -  x2 
x=l  m 

4  /  /  x1  y]  dy  dx  , 

x=0  y=0 

where  the  proper  i  and  j  are  taken  from  the  term  being  integrated.  Integration 
with  respect  to  y  simplifies  this  expression  to: 


4 

j+1 

(j  +  l)nv 


o 


The  evaluation  of  Q0  and  Q2  is  reduced  to  the  evaluation  of  this  expression  for 
each  term  in  the  equations,  multiplying  it  by  the  corresponding  constant  and 


summing  the  values  obtained  for  each  term  in  the  equation. 

A  computer  program  was  written  for  the  preceding  W-,  X-,  Y  ,  and 
Z-matrices  and  the  doable  integration  ot  the  eolations.  The  following  20  cases 


are  presented: 


Case  1 


n  =  0. 1 


C  =  2.02 


~f~=  1  -  0  .  785  935  X  10“9K2 
s 


Case  2 


_  ?rWna2 


fl  -  0.785935  X  10"9K2] 


m  =  0.2 


C  *  2.08 


— =  J  "  0.140091  X  10“7K2 
s 


Case  3 


QT  "  -  Q-J-  f1  ~  0.140091  X  10"7K2] 


m  =  0.3 


C  =  2.18 


— =  1  -  0.623529  X  10_7K5 


C  se  4 


qt  =  .1^  fi  .  o. 


623529  X  10~7k2J 


m  =  0. 4 


C  -  2.32 


P  -  1  -  0.  136545  X  10~6k2 
s 


7rW„a 


-  M  -  0.  135545  X  lO-^K2] 
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Case  5  m  =  0.5  C  =  2.50 

F 

r~  «  1  -  0. 194319  X  10-6K2 
s 

Qt  =  7rW0a2 [  1  -  0. 19431b  X  10"eK2] 

Case  6  m  =  0.6  C  =  2.72 

F 

-^-  =  1  -  0.209387  X  10“6K2 
s 

_  TrW'^a2  _  209387  x  10-e  :2j 

T  i.2 

Case  7  m  =  0.7  C  =  2.98 

F 

■jr  =  1  -  0. 188178  X  10-6K2 
s 

q  =  7r^..°a  .  [i  .  o.  188178  X  10"eK2 
T  1.4 

Case  8  m  =  C.8  C  =  3.28 

F 

=  1  -  0. 150800  X  10~eK2 
s 

Q  [1  -  0.150800  X  10~eK2] 

T  1.6 

Case  9  m  =  0.9  C  =  3.62 

F 

-“=  1  -  0. 112857  X  10-eK2 
*  s 

Q  =  ?Wq...-  [1  -  0. 112857  X  10-GK2j 

f  1.8 
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Case  15 

F 

m  =  1.5 

C  =  6,50 

> 

: 

/  =  1 

-  0. 143396  X  10"7K2 

s 

' 

Q 

T 

_  7rW0a2 

3 

[1  -  0.143396  X  10-7K2] 

Case  16 

F 

m  -  1.  6 

C  =  7.12 

T^1 

-  0.102846  X  10-7K2 

S 

| 

Q 

T 

_  7rW0a2 

[1  -  0.102846  X  10-7K2] 

[ 

3.2 

Case  17 

F 

m  =  1.7 

C  =  7.78 

/-1 

s 

-  0. 743664  X  10“8K2 

Qm 

T 

7rW0a2 

[1  -  0.743684  X  10~8K2] 

3.4 

Case  18 

F 

m  =  1,8 

C  =  8.48 

1 

-  0.542271  X  10"8K2 

!  ‘ 

s 

$ 

Qrp 

7rW0a2 

[1  -  0.542271  X  10"8K2] 

3.6 

Case  19 

F 

m  =  1. 9 

C  =  9.22 

j 

T"=  1 

-  0.398714  X  10“8K2 

1 

s 

j 

qt 

7rW0a2 

3.8 

[1  -  0.398714  X  10"8K2] 

; 

|  • 

i  : 

j 

i 
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Case  20 


m  =  2.0 


C  =  10.00 


F 

—  -  1  -  0.295576  X  10"8K2 
s 


rW0a2 

4 


[1  -  0.295576  X  10“SK2] 


■  tituting  equation  (4.6)  into  (4.1)  and  integrating  yields: 


Qo  _ 


7rW0a?C 
4m(m2  +  1) 


(4.8) 


which  is  the  rate  of  flow  through  a  straight  pipe  with  an  elliptical  cross  section 
for  any  C  >  0  and  m  >  0. 

The  computer  program  solves  the  problem  for  any  C  >  0,  but  for  the 
examples  presented  C  =  2(m2  +  1)  by  choice.  If  m  =  1;  c  =  4,  which  is  the 
value  of  C  taken  by  Dean  for  his  solution  and  equation  (4.8)  reduces  to: 


2m 


(4.9) 
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CHAPTER  V 

DISCUSSION 

After  formulation  of  the  present  theory,  a  digital  computer  program  was 
written  to  obtain  the  required  solutions.  The  program  is  in  Appendix  A .  To 
get  the  rate  of  flow  through  a  pipe  with  a  small  curvature  and  an  elliptical 
cross  section  requires  the  following  input  data:  (1)  the  value  of  C  for  the  Puid, 
(2)  the  value  of  m  for  the  elliptical  cross  section,  and  (3)  the  value  of  K 
(Dean's  number) . 

Twenty  cases  were  considered  for  C  =  2(m2  +  1) ,  m  =  0. 1,  0.2, 

0. 3,  ....  2.0.  The  solutions  were  expressed  in  terms  of  K  and  the  nondimen- 
sionalization  constants  W0  and  a.  Tne  results  for  these  cases  are  presented  in 
Chapter  IV  and  Appendix  B. 

The  graphs  in  this  section  are  plotted  for  the  cases  C  =  2(m2  +  1)  and 
m  =  0.5,  1.0,  and  1.5.  The  cross  section  of  a  pipe  is  shown  in  Figure  3  for 
each  value  of  m. 

A .  Streamlines  in  the  Cross  Sectional  Plane  of  the  Pipe 
and  the  Vorticity  Centers  of  the  Secondary  Flow 

The  first-order  approximation  of  the  streamlines,  ^  =  constant,  was 
plotted  by  Thomas  and  Walters  [1964]  for  a  pipe  with  an  elliptical  cross  section. 
The  second-order  approximation  of  the  streamlines,  ip2  =  constant,  has  not 
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t 

previously  been  plotted  for  a  pipe  with  an  elliptical  cross  section.  The  centers 
of  secondary  flow  were  presented  by  Dean  [1927]  for  a  pipe  with  a  circular 
cross  section.  The  following  is  an  extension  of  the  existing  data  and  is  in 
complete  agreement  with  it. 

1.  First-Order  Approximation  of  Streamlines 
on  the  Cross-Sectional  Plane  of  the  Pipe 

To  plot  the  streamlines,  =  constant,  on  a  cross  section  of  the  pipe 

equation  (3.2)  is  expressed  as 

[Ajyjx6  +  [  (A0  -  2Aj)y  +  (2Ajm2  +  A2)y3]x4  (5.1) 

+  [(A,  -  2A0)y  +  ( 2A0m2  -  2A,m2  -  2A2)y3 
+  (Ajtn4  +  2A2m2)y5]x2  [A0y  +  (A2  -  2A0m2)y3 
+  (A0m4  -  2A2m2)y5  +  A2m4y7  -  <//]  =  0  . 

Constant  values  are  taken  for  ipi  and  y.  This  reduces  the  bracketed  terms  to 
constants  resulting  in  a  sixth-degree  polynomial  in  x  with  constant  coefficients. 
The  data  required  for  a  plot  of  x  versus  y  with  i ^  =  constant  are  calculated  by 
subroutine  PLOT  of  the  digital  computer  program  in  Appendix  A.  These  data 
are  presented  in  Tables  1,  2,  and  3  and  plotted  in  Figures  4,  5,  and  6. 

Figures  4,5,  and  6  show’,  as  can  be  seen  from  equation  (5.1),  that  the  value 
of  x  is  undetermined  if  y  =  0.  The  streamlines,  i ^  =  constant,  are  symmetric 
with  respect  to  the  /-axis.  Taking  the  negative  value  of  the  constant  used  to 
plot  the  streamlines  in  the  first  and  second  quadrants  results  in  another  curve 
in  the  third  and  fourth  quadrants.  The  streamlines  plotted  for  the  two  constants 


swss&ai 
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Table  1.  First-Order  Approximation  of  Streamlines 
in  the  Cross  Section  of  the  Pipe  at  m  =  0.5 


y 

X 

|,  =  5X  10“4 

X 

’l‘i  =  2  X  10“3 

! 

\ 

1.7 

0.2037 

-0.2037 

! 

1 

1.3 

0.5673 

-0.5673 

0.2966  -0.2966 

£ 

■  i 

j 

0.9 

0.7129 

-0.7129 

0.4902  -0.4902 

* 

i 

$ 

0.5 

0.7551 

-0.7551 

0.4674  -0.4674 

t 

i 

0.1 

0.4538 

-0.4538 

j 

ipt  =  -5  X  10' 


-0.4538 

-0.7551 

-0.7129 

-0.5673 

-0.2037 


0.4538 

0.7551 

0.7129' 

0.5673 

0.2037 


-0.4674  0.4674 

-0.4902  0.4902 

-0.2966  0.2966 


are  symmetric  about  the  x-axis  and  the  origin.  The  x,y  relation  is  dependent 
on  the  parameters  C  and  m.  In  the  cases  plotted  C  =  2(m2  +  1)  by  choice.  In 
loose  terms,  the  streamlines  may  be  thought  of  as  projections  of  the  paths  of 
fluid  elements  on  the  cross  section  of  the  pipe.  The  secondary  motion  is  caused 
by  the  pipe  curvature  and  represents  a  loss  of  energy  which  retards  the  primary 
flow  through  the  pipe  for  a  given  inlet  pressure. 


Table  2.  First-Order  Approximation  of  Streamlines 
in  the  Cross  Section  of  the  pipe  at  m  =  1.0 


01  = 

X 

1  X 

10~4 

= 

X 

7  X 

10"4 

0.2175 

-0.2175 

0.5905 

-0.5905 

0.3147 

-0.3147 

0. 7483 

-0.7483 

0.5184 

-0.5184 

0.8165 

-0.8165 

0.5475 

-0.5475 

0.7610 

-0.7610 

= 

-1  X 

10"4 

01  = 

-7  X 

io-4 

-0.7610 

0.7610 

-0. 8165 

0.8165 

-0.5475 

0.5475 

-0.7483 

0.7483 

-0.5184 

0.  5184 

-0.5905 

0.5905 

-0.3147 

0.3147 

-0.2175 

0.2175 

Second-Order  Approximation  of  Streamlines 


on  the  Cross-Sectional  Plane  of  the  Pine 


To  plot  the  streamlines  (ip2  =  constant)  on  a  cioss  section  of  pipe, 
equation  (3.4)  is  expressed  as 


l-ioy]xi:)  +  l(-2C10  +  C6)y  +  (C„  +  2m2C10)y3]x11  +  [  (C10  -  2CG  +  C3)y 

(5. 

+  (~2Cu  -  2C]0m2  +  C7  ^  2m2C6)y:i  +  (C]2  +  2Cj<m2 


+  m4C10)y5lx»  +  [(C8  -  2C3  *  C,)y  (C1t  -  2C7  -  2C6m2  +  C4 


Table  3.  First-Order  Approximation  of  Streamlines 
in  the  Cross  Section  of  the  Pipe  at  m  --  1.5 


y 

X 

^,  =  5X 

10“5 

, 

01  = 

X 

3  X 

*0* 

l 

O 

1 

0.5 

0.1648 

-0.1648 

0.4 

0.2394 

-0.2394 

0.4171 

-0.4171 

0.3 

0.3538 

-0.3538 

0.5073 

-0.5073 

0.2 

0.3195 

-0.3195 

0.5077 

-0.5077 

0. 1 

0.3158 

-0.3158 

0j  -  -5  X 

o 

1 

0i  =  ■ 

-3  X 

V— * 

o 

1 

0 

-0.1 

-0.3158 

0.3158 

-0.2 

-0.3195 

0.3195 

-0.5077 

0.5077 

-0.3 

-0.3538 

0.3538 

-0.5073 

0.5073 

-0.4 

-0.2394 

0, 2394 

-0.4171 

0.  4171 

-0.  5 
. . . .  _ 

-0. 1648 

0.1648 

+  2irTC:))y,i 

<  (-2Cj2  -  2C1Jnr  -i 

C8  +  2nrC7  +  mCe)  y 

+  2nrC12  ‘ 

ni'1C1i)y7ix7  f(C3  - 

2C,  +  C0)y  +  (C7  - 

-  2C3nv’  • 

C,  ‘  2m2C,)y3  +  (C-2 

-  2Cg  -  2C7m2  +  C 

*  m'Cjly'1 

■  (-2C13  -  2C,,m2  t  ( 

|  2m?C8  1  m4C7)y7 

+  2ni2C13  i 

m4Cl2)y‘'!x5  >  [  ( C ,  - 

2C0)y  *  (C,  -  2 C,  - 

2inzC0)y3  '  (C„  -  2C5  -  2C.,nr  *  2m2C,  -  m4C,)v'’ 


i 
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+  m4C8)y9  +  (2m2C14  +  m'c^y^Tx3  +  [Ctf  +  {C2  -  2Com2)y3 
+  (Cs  -  2C2m2  +  m4C0)y5  +  (C9  -  2Csm2  +  C2m4)y7  +  (C14 
-  2C9m2  +  m4C5)y9  +  (-2C14m2  +  m4C9)y11  +  m4C14y13Jx  -  tp2  =  0 

Constant  values  are  taken  for  ip2  and  y.  This  reduces  the  bracketed  terms  to 
constants  resulting  in  a  thirteenth-degree  polynomial  in  x  with  constant  coeffi¬ 
cients.  The  data  required  for  a  plot  of  x  versus  y  with  $2  =  constant  are  cal¬ 
culated  by  subroutine  PLOT  of  the  digital  computer  program  in  Appendix  A. 
These  data  are  presented  in  Tables  4,  5,  and  6  and  plotted  in  Figures  7,  8,  and 
9.  Figures  7,  8,  and  9  show,  as  can  be  seen  from  equation  (5.2) ,  that  the 
value  of  x  is  undetermined  if  y  =  0  and  the  value  of  y  is  undetermined  if  x  =  0. 
The  curves,  tp2  =  constant,  are  plotted  in  the  first  and  third  quadrants  by  taking 
a  positive  constaut  for  </i2,  and  taking  a  negative  value  of  the  same  constant  pro¬ 
duced  the  same  curves  in  the  second  and  fourth  quadrants.  The  curves  in  the 
four  quandrants  are  symmetric  about  the  x-axis,  y-axis,  and  the  origin.  The 
x,y  relation  is  dependent  on  the  parameters  C  and  m.  In  the  cases  plotted 
C  =  2(m2  +  1)  by  choice.  It  is  clear  from  Figures  7,  8,  and  9  that  the  contri¬ 
bution  of  $2  to  the  secondary  velocities,  U  and  V,  are  the  same  for  each  of  the 
four  quadrants. 

3.  Vorticity  Centers  of  the  Secondary  Flow 

The  centers  of  secondary  flow  are  points  in  the  cross-sectional  plane  of 
the  pipe  where  the  secondary  velocity  vanishes.  These  points  are  called 
vorticity  centers.  For  the  secondary  velocity  to  vanish  both  components, 


Table  4.  Second- Order  Approximation  of  Streamlines 
in  the  Cross  Section  of  the  Pipe  at  m  =  0.5 


h  = 

X 

3  X 

io-7 

h  = 

X 

5  X 

io-8 

- 

0.2427 

0.1902 

0.5709 

0.0284 

0. 7883 

0.1275 

0.7250 

0.205 

0.8837 

0.4286 

0.4286 

0.0740 

0. 8243 

-0.4286 

-0.4286 

-0.0740 

-0. 8243 

-0.1276 

-0.7250 

-0. 0205 

-0.8837 

-0. 1902 

-0.5709 

-0.0284 

-0.7883 

-0.2427 

h  = 

-3  X 

io-7 

h  = 

-5  X 

1 

CO 

o 

rH 

-0.2427 

-0.1902 

-0.5709 

-0.0284 

-0.7883 

-0.1276 

-0.7250 

-0.0205 

-0.8837 

-0.4286 

-0.4286 

-0.0740 

-0.8243 

0.4286 

0.4286 

0. 0740 

0.8243 

0.1276 

0.7250 

0.  0205 

0.8837 

0.1902 

0.5709 

0.0284 

0.7883 

0.2427 

Table  5.  Second-Order  Approximation  of  Streamlines 
in  the  Cross  Section  of  the  Pipe  at  m  =  1.0 


lp2  - 

X 

5  X 

io-8 

h  - 

X 

1  X 

io-7 

0.1998 

' 

0.5778 

0.5778 

0.1519 

0.4633 

0.0387 

0.7409 

0. 0787 

0.6787 

0.0400 

0.  8164 

0.0812 

0.7492 

0.0989 

0.7576 

0. 6070 

-0.0989 

-0.7576 

-0. 6070 

-0.0400 

-0. 8164 

-0.0812 

-0. 7492 

-0.0387 

-0. 7409 

-0.0787 

-0. 6787 

-0.5778 

-0,5778 

-0.1519 

-0.4633 

-0.1998 

h  = 

-5  X 

io-8 

h  = 

-1  X 

10"7 

-0.1998 

-0.5778 

-0.5778 

-0.1519 

-0.4633 

-0.0387 

-0.7409 

-0. 0787 

-0.6787 

-0.0400 

-0.8164 

-0.0812 

-0.7492 

-0.0989 

-0.7576 

-0. 6070 

0,0989 

0.7576 

0.6070 

0.0400 

0.8164 

0. 0812 

0. 7492 

0  .  03  87 

0.7409 

0.0787 

0.6787 

0.5778 

0.5778 

0.1519 

0.4633 

0.  1998 

Table  6.  Second-Order  Approximation  of  Streamlines 
in  the  Crosfs  Section  of  the  Pipe  at  m  =  1.5 


II 

CM 

_ 1 

X 

5  X 

10~3 

h  = 

X 

2  X 

10"8 

0.2703 

0.1056 

0.4569 

0.1927 

0.4585 

0.0678 

0.6180 

0. 1644 

0.5572 

0.0517 

0.6921 

0.1705 

0.5902 

0.0356 

0.7254 

0.4611 

0.4786 

-0.4611 

-0.4786 

-0.1705 

-0.5902 

-0.0356 

-0.7254 

-0. 1644 

-0.5572 

-0.0517 

-0.6921 

-0.1927 

-0.4585 

-0.0678 

-0.6180 

-0.2703 

-0.1056 

-0.4569 

^2  =  • 

-5  X 

10"8 

h  =  • 

-2  X 

H- ■* 

o 

1 

CO 

-0.2703 

-0.1056 

-0.4569 

-0. ‘ 927 

-0.4585 

-0.0678 

-0. 6180 

-0.1644 

-0.5572 

-0.0517 

-0.6921 

-0.1705 

-0.5902 

-0.0356 

-0. 7254 

-0.4611 

-0.4786 

0.4611 

0.4786 

0.1705 

0.5902 

0. 0356 

0.7254 

0.1644 

0.5572 

0.0517 

0.6921 

0.1927 

0.4585 

0.0678 

0.6180 

0.2703 

0.1056 

0.4569 
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U  and  V,  must  vanish.  Sufficient  accuracy  for  location  of  the  vorticity  centers 
is  retained  after  the  second-order  approximation  of  the  stream  function  is 
neglected.  Applying  equations  (2.34)  and  (2.37),  and  differentiating  equation 
(3.2) ,  yields: 

V  =  (1  -  x2  -  m2y2)2xy[  (2A0  -  Aj)  +  3AjX2  +  (Ajm2  +  2A2)y2] 

dx  (5.3) 

U  =  =  (1  -  x2  -  m2y2)  [A0  +  (Aj  -  A0)x2  +  (3A2  -  5A0m2y2  (5.4) 

3y 

+  (-Aj)x4  +  (-SAjm2  -  3A2) x2y2  +  (-7A2m2)y4] 

Inspection  of  equation  (5.3)  reveals  that  V  vanishes  for  all  values  of  y  at  x=  0. 
Setting  U  =  0  in  equation  (5.4)  dividing  by  the  boundary  conditions  and  substi¬ 
tuting  x  =  0  yields: 

7A2m2y4  +  (5A0m2  -  3A2)y2  -  A0  =  0  .  (5.5) 

Solving  equation  (5.5)  by  the  quadratic  formula  and  taking  the  square  root  of 
both  sides  yields: 


The  minus  sign  on  the  small  radical  is  removed  because  it  results  only  in 
imaginary  values  of  y,  which  are  not  applicable  to  this  problem.  The  data 
required  for  a  plot  of  m  versus  y  are  calculated  by  the  digital  computer  pro¬ 
gram  in  Appendix  A.  These  data  are  presented  in  Table  7  and  plotted  in 


> 


Table  7.  Position  of  the  Vorticity  Centers  of  the  Secondary  Flew 


Figure  10.  The  curves  in  Figure  10  are  symmetric  about  the  m-axis.  The 
m,y  relation  is  dependent  on  the  parameters  C  and  m.  In  the  cases  plotted 
C  =  2(m2  +  1}  by  choice.  It  is  evident  from  Figure  10  that  y  approaches  zero 
as  m  becomes  infinite  and  y  approaches  infinity  as  :m  approaches  zero.  This  is 
to  be  expected  because  the  axis  of  the  ellipse  coinciding  with  the  y-axis,  B,  is 


STw¥i 


[I] 
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determined  by  the  relation  B  =  l/m  and  the  y  coordinate  of  the  vorticity  centers 
increases  with  an  increase  in  B.  For  m  =  1  (circular  case)  the  position  of 
the  vorticity  center  is  the  same  as  that  obtained  by  Dean  [1928]. 


B.  Streamlines  in  the  Central  Plane 
The  dimensional  differential  equation  for  the  streamlines  in  the  central 
plane  is: 


(R  =  x')d0 


q 
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dxT 


(5.7) 


Sufficient  accuracy  is  retained  if  only  the  zero-order  approximation  of  W  and 
the  first-order  approximation  of d  are  considered.  Also,  xr  is  negligible  in 
comparison  with  R  for  small  curvatures.  At  the  central  plane  y  =  0,  and  from 
equations  (2.24),  u,31),  and  (2.33)  the  terms  become: 


rU  Svj  rK  2  2  i 

G  =  — =  —  — —  = -  fl  -  xz)z(A0  ~  AiXT)  at  v  =  0 

a  a  ay  a  q  1 


(5.8) 


dx'  =  adx 


q  =  W0(l  -  x  )  at  y  =  0 

v 


where  C  =  2(m2  4-  1)  by  choice.  Substituting  equations  (5.8)  into  equation 
(5. 7)  and  solving  for  de  in  nondimensional  terms  yields: 


r 


de  = 


2W0a 


dx 


(1  -  x2)  (A0  +  AjX2) 


(5.9) 
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Integration,  of  equation  (5.9)  yields: 


*  *  / — I —  1  +  x  /  — — 

- i -  logi±JU  /-illog _ N  A» 

— (Ao  +  Ai>  1  -  x 

v  Aa 


(5.10) 


when—  <  0. 
A0 


1  ,  1 +  x  ■  /Ai  t  /  /Ai 

4^(a0  +  Ai)  ^  W'A« 

r  L 


(5. 11) 


A, 

when — i  >  0. 
A0 


The  data  required  for  a  polar  plot  of  9  versus  x  are  calculated  by  sub¬ 
routine  PLOT  of  the  digital  computer  program  in  Appendix  A .  The  term  9 
increases  steadily  with  x  and  approaches  infinity  as  x  approaches  one  or  a  mini  s 

one.  Therefore,  the  streamline  approaches  but  never  reaches  the  inside  sur- 

a 

face  of  the  pipe.  The  relation  between  9  and  x  does  not  involve  —  and  is.  there- 

R 

fore,  independent  of  the  curvature  of  the  pipe.  Since  A0  and  depend  on  the 
parameters  C  and  m ,  the  variation  of  9  with  x  depends  on  C ,  m,  and  the 
Reynolds  number.  In  Table  S  the  calculated  values  of  9  versus  x  are  given  for 
in  =  0.5,  1.0,  and  1.5,  taking  the  value  of  4W0a/r  to  be  130.  The  form  of 
streamlines  for  ordinary  flow  are  shown  in  Figure  11.  In  Figure  11  —has  been 

a 

assumed,  for  graphical  illustration,  to  have  a  large  yalue  of-|  since  the  x,  0 

g 

relation  is  independent  of—,  but  the  data  would  not  be  valid  for  a  pipe  with  such 

K 


large  curvature.  From  Figure  11  it  can  be  seen  that  the  curvature  of  the 


Table  8.  Streamlines  in  the  Central  Plane 


X 

0 

m  =0.5 

6 

m  =  1.0 

e 

m  =  1.5 

1.0 

or 

oo 

oo 

0.8 

69.1 

110.9 

312.0 

0.6 

44.1 

71.7 

204.0 

0.4 

27.1 

44.4 

126.2 

0.2 

13.0 

21.4 

60.9 

0 

0 

0 

0 

-0.2 

-13.0 

-  21.4 

-  60.9 

-0.4 

-27.1 

-  44.4 

-126.2 

-0.6 

-44.  J 

-  71.7 

-204.0 

-0.8 

-69.1 

-110.9 

-312.0 

•1.0 

_ 

OC 

oc 

OO 

streamline  in  the  central  plane  decreases  as  m  increases. 

C.  Effe.  ;  of  Pipe  Curvature  on  the  Flow  Rate 
The  effect  of  pipe  curvature  on  the  flow  rate  is  given  by  equation  (4.5) . 


It  is 


F  r 


s  L_ 


1  -  2*  K: 


Qo 


where  F  /F  is  the  ratio  of  the  flux  through  a  curved  pipe  to  the  flux  through  a 

C  3 

straight  pipe,  bothhavingthe  same  inlet  pressure,  length,  and  cross  section.  The 
Q? 

ratios,  _  ,  are  calculated  for  20  values  of  m  by  the  digital  computer  program 

Qo 

in  Appendix  A.  The  ratios,  F  /F  ,  are  calculated  for  values  of  K  ranging  from 

c  s 


0  to  1000  by  jse  of  a  desk  calculator.  These  data  are  presented  in  Table  9  and 
plotted  In  Figure  12  for  m  =  0. 5, 1. 0,  and  1. 5.  The  F J Fg,  K  relation  is  depen¬ 
dent  on  the  parameters  C  and  m  since  Q2  and  Qc  are  dependent  on  C  and  m.  In 

the  cases  plotted  C  =  2(m2+  1)  by  choice.  The  F  /F  values  taken  from 

c  s 

Figure  12  ire  corrections  factors  for  the  pipe  curvature.  The  F  /F  values 

c  s 

are  used  with  the  equation 


(5.12) 


where  Q  is  the  total  flow  rate  through  a  curved  pipe  and  Qq  is  the  flow  rate 

through  i  straight  pipe  with  the  same  inlet  pressure,  length,  and  cross  section. 

From  Figure  12  it  is  noted  that  F  /F  is  approximately  one  for  an  extremely 

c  s 

small  Dean's  number  and  the  curvature  effect  is  negligible.  As  the  Dean's 
number  i:  increased,  the  curvature  effect  rapidly  increases.  This  is  expected 
since  the  Dean's  number  is  proportional  to  the  velocity  squared  and  the  energy 
loss  is  much  greater.  The  curvature  effect  decreases  with  an  increase  in  the 
value  of  m.  Figure  12  also  shows  that  the  effect  of  pipe  curvature  on  the  rate 
of  flow  through  the  pipe  is  greater  for  a  pipe  with  the  major  axis  of  the  cross- 
sectional  ellipse  perpendicular  to  the  plane  of  the  bend  (m  <  1)  than  when  the 
major  axi.i  coincides  with  the  radius  of  curvature  (m  >  1) . 

D.  First-Order  Approximation  of  the  Primary  Velocity 
in  the  Central  Plane  ~~ 

The  first-order  approximation  of  the  primary  velocity  is  given  by 
equation  (3.3)  as 


J 
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Table  9.  Flux  Through  a  Curved  Pipe/Flux  Through  a  Straight  Pipe 
for  Various  Dean's  Numbers 


K 

F  /F 
c  s 

m  =  0.5 

F  /F 
c  s 

m  =  1. 0 

F  /F 
c  s 

m  =  1.5 

0 

1.000 

1.000 

1.000 

10 

0.999 

0.999 

0.999 

100 

0.998 

0.999 

0.999 

200 

0. 992 

0.997 

0.999 

300 

0.983 

400 

0.968 

0.987 

0.997 

500 

0.951 

600 

0.930 

0.971 

0.995 

700 

0.905 

800 

0.876 

0.948 

0.991 

900 

0.843 

1000 

0.806 

0.919 

0.986 

Wj  =  (I  -  x2  -  m2y2)x[b0  +  bjx2  +  b2y2  +  b3x4  +  b4x2y2 
+  bjy4  +  b6x6  +  b7x4y2  +  b8xV  +  b9y6] 


In  the  central  plane  y  =  0  and  equation  (3.3)  becomes: 


Wj  =  b()X  +  (bj  -  b0)x3  +  (b3  -  b^x5  +  (b6  -  b3)x7  -  b6x9  (5.13) 

The  data  required  for  a  plot  of  x  versus  Wj  are  calculated  by  subroutine  PLOT 
of  the  digital  computer  program  in  Appendix  A.  The  factor  104  is  just  an 
amplification  factor  for  the  curves.  These  data  are  presented  in  Table  10  and 
plotted  for  m  =  0.5,  1.0,  and  1.5  in  Figure  13.  Figure  13  shows,  as  can  be 


400  600 

DEAN’S  NUMBER,  K 


800 


1000 


of  Pipe  Curvature  Versus  Dean's  Number 
th  m  =  0. 5,  1.0,  1.5 


Table  10.  First  Order  of  Approximation 
of  the  Primary  Velocity  at  y  =  0 


X 

Wj  X  104 

W,  X  1G4 

W,  X  104 

- 

m  -  0.5 

m  =  1.  0 

1.0 

0 

0 

I 

0 

0.8 

!  1.897 

1.123 

0.431 

0.6 

3.336 

2.093 

0.877 

0.4 

3.546 

2.314 

1.026 

0.2 

2.271 

1.515 

0.693 

0 

0 

0 

0 

-0.2 

-2.271 

-1.515 

-0. 693 

-0.  4 

-3.546 

-2.314 

-1.026 

-0.6 

-3.336 

-2.093 

-0. 877 

-0.3 

-I. 897 

-1.123 

-0.431 

-1.0 

0 

o 

_ i 

0 

seen  from  equation  (5. 13) ,  that  the  function  is  odd  and  symmetric  about  the 
origin.  The  x,  Wj  relation  is  dependent  on  the  parameters  C  and  m.  In  the 
cases  plotted  C  =  2(m2  +  1)  by  choice.  Figure  13  shows  that  the  magnitude  of 
the  f^rst -order  approximation  of  the  primary  velocity  in  the  central  plane 
decreases  with  an  increase  in  the  value  of  m. 
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CHAPTER  VI 

SUMMARY  AND  CONCLUSIONS 

The  graphical  illustrations  presented  show  that  the  streamline  flow  of  ' 
an  incompressible  fluid  through  a  curved  pipe  of  elliptical  cross  section  is 
similar  to  that  of  a  curved  pipe  of  circular  cross  section.  It  consists  of  a 
primary  flow  and  a  secondary  flow.  The  primary  flow  is  along  and  parallel  to 
the  center  line  of  the  pipe.  The  secondary  flow  is  in  the  plane  of  the  cross 
section  of  the  pipe  with  the  form  of  the  streamlines  as  shown  in  Figures  7 
through  9.  The  motion  of  the  secondary  flow  ineffectively  expends  energy  .which 
results  in  a  decreased  rate  of  flow  through  the  pipe.  The  secondarjTflowand 
centrifugal  forces  acting  on  the  fluid  are  not  present  in  a  straight  pipe.  There¬ 
fore,  the  rate  of  flow  through  a  straight  pipe  is  an  upper  bound  for  the  rate  of 
flow  through  a  curved  pipe  with  the  same  inlet  pressure,  length  and  cross 
section. 

The  assumption  made  in  the  derivation  of  the  governing  differential 
equations  are  equivalent  to  those  made  by  Dean  [1928]  for  a  pipe  with  a  circular 
cross  section.  The  governing  differential  equations  of  this  study  are  equivalent 
to  those  obtained  by  Dean  [1928]  and  the  accuracy  is  the  same.  However,  the 
solutions  to  the  equations  obtained  in  this  study  are  exact  and  applicable  to  a 
pipe  with  an  elliptical  cross  section.  Dean's  [1928]  solutions  are  approximate 
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because  the  unsymmetrical  terms  were  neglected  in  the  equation  for  the  second- 
order  approximation  of  the  primary  velocity. 

For  the  circular  case  (m  =  1)  the  solutions  to  the  first  four  of  the 
governing  differential  equations  (2.38)  through  (2.41)  are  exactly  equivalent 
to  the  corresponding  solutions  presented  by  Dean  f  1928]  and  are  symmetric. 

The  solution  to  the  fifth  governing  equation  (2.42)  differs  from  Dean's  [1928] 
approximate  solution  by  the  value  represented  by  the  unsymmetrical  terms. 

The  error  in  Dean's  [1928]  second-order  approximation  of  the  flow  rate,  due 
to  neglect  of  the  unsymmetrical  terms,  is  approximately  0. 17  per  cent  and 
1.2  per  cent  for  K  values  of  400  and  1000,  respectively.  The  flow  rate  given  by 
the  equations  of  the  present  study  differs  from  that  of  Dean's  [1928]  fourth- 
order  approximation  by  0.12  per  cent  and  0.  0044  per  cent  for  K  values  of  400 
and  1000,  respectively.  Therefore,  the  equations  obtained  in  the  present  study 
have  about  the  same  accuracy  as  Dean's  [1928]  fourth-order  approximation  for 
the  case  of  a  circular  cross  section. 

Figure  12  shows  that  the  effect  of  pipe  curvature  on  the  rate  of  flow  is 
greater  for  a  pipe  with  the  major  axis  of  the  cross-see. ional  ellipse  perpendic¬ 
ular  to  the  plane  of  the  bend  (m  <  1)  than  when  the  major  axis  coincides  with 
the  radius  of  curvature  (m  >  1) .  However,  the  fabrication  is  generally  more 
difficult  for  pipe  bends  with  m  >  1. 

The  rate  of  flow  of  fluid  through  a  curved  pipe  with  an  elliptical  cross 
section  is  a  function  of  C,  m,  and  K.  Therefore,  the  rate  of  flow  of  the  fluid 
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ia  dependent  on  the  fluid  properties ,  cross  section  of  the  pipe  bend  and  the 
Dean's  number. 

The  pipe  curvature  reduces  the  rate  of  flow  through  the  pipe  if  the  inlet 
pressure,  pipe  length,  radius  of  curvature,  end  cross  section  of  the  pipe  are 
identical.  The  magnitude  of  this  reduction  decreases  as  the  value  of  m 
increases  (Figure  12).  Therefore,  the  accuracy  of  the  solutions  presented  is 
dependent  upon  m,  but  the  governing  equations  are  probably  net  valid  for  any 
m  when  the  Dean's  -lumber  exceeds  1000. 


) 
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CHAPTER  VH 

RECOMMENDATIONS  FOR  FUTURE  RESEARCH 

The  following  recommendations  are  made  for  fixture  research: 

1.  A  more  accurate  solution  can  probably  be  obtained  by  transferring 
the  governing  equations  (2.35)  and  (2.36)  into  finite  differei  ne  form  and 
solving’  them  directly  on  a  computer. 

2.  The  accuracy  of  this  solution  can  be  impro/ed  by  extending  it  to  a 
fourth-order  approximation,  but  considerable  work  will  be  involved. 

3.  The  solutions  presented  could  be  generalized  to  cover  clastico- 
viscous  fluid  or  magnetobydrodynamic  flow. 

4.  The  solutions  presented  could  be  further  verified  by 
experimentation . 
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APPENDIX  A 

OPERATING  INSTRUCTIONS 
AND 

COMPUTER  PROGRAM 

This  program  is  for  use  on  the  IBM  7094  Computer.  To  solve  for  the 
rate  of  flow  through  the  pipe  and  generate  the  data  required  to  plot  the  graphs 
presented  in  Chapter  V,  numerical  values  for  the  case  cumber,  XM,  C,  XK, 

PSI  11,  PSI  21,  PSI  12,  PSI  22,  PSI  13,  PSI  23,  PSI  14,  PSI  24,  Y0,  and  DEL  Y 
must  be  included  as  input  data.  The  input  data  are  defined  as  follows: 

XM  =m 

C  =  C 

XK  =  K  =  Dean's  number 

PSI  11,  PSI  12,  PSI  13,  and  PSI  14  -  four  constants  taken  for 

PSI  21,  PSI  22,  PSI  23,  and  PSI  24  -  four  constants  taken  for  #2 

Y0  =  initial  value  of  Y  for  which  values  of  X  are  determined 

DEL  Y  =  increment  of  Y  added  to  Y3 

If  the  data  to  plot  the  graphs  are  not  desired,  the  CALL  PLOT  card  with 
its  two  continuation  cards  at  call  number  1004  can  be  removed  from  the  deck 
and  the  program  ends  after  solving  for  the  rate  of  flow  through  the  pipe.  How¬ 
ever,  values  for  all  of  the  input  variables  must  be  included  as  input  data.  Since 

the  only  input  variables  used  to  determine  the  flow  rate  are  XM,  C,  and  XK, 
any  constant  can  be  uQed  for  the  other  input  data. 


IEJC9  VERSICA  5  HAS  CCATFOL. 
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AIEFTC  PAIA  LIST. REF 


COUBLfc  PRECISION  M.X.Y.Z.OET 

DOUBLE  PRECISION  XK.XM2, AO,  A  i.A2, 80,  Hi, 6  2. 83, 8 .,85 ,B6 ,87 , B8 ,89 ,CO  , 
lCl.C2.C3,CA,C5,Cfc,C7,C8,CS,ClQ.Cll,C12.Cl3,ClA,O0,01,D2,G3,  OA.05,0 
26. 07 . OE. 09. C 10. Oil, D12, 0 1?, Cl A, 015, 016, 01 7,016, 019, D20. 021, 022, 023 
3.C2A  .025 ,026 ,027. 028 .029,1)  30,031  ,032,031 ,03 A  ,D35,C,A 
0 1  HE  NS  1C  N  M 1 36 ,37 1 
CIREASICN  X<?6,37> 

CINE  AS  ICN  Yl  36,37) 

DINEASICA  2(36,27) 

CCRNCH  XI.XJJ 

eoui valeace  iwu.ii.xn.n.Yi l.n.zi l.i) i 
INTECER  CASE 

NAME  L 1  ST  /  INP1  T  /CASE  ,  XP  ,C  PS  I  U ,  PSI  2 1  ,PS  I  12 ,  RSI  22 ,  PS  1 1  3,  PS!  23  ,PS  1 1 A 
1.PS!2a,YC.XK,0EIY 
S  PEACIS. INPUT) 

KRiTEIC.  1APLT1 

C  =  C.C 

XH2  *  XJ**XM 

A=C/  I2.*IXM**2*1  .11 

tot  1,  1)  =  5.*XH**A*2 .*XH**2U  . 

*! 1.2  1 =-2 .*XP**2-2. 

W* 1.  2 ) *- 1C. *XN**2- 2 . 

h 1 l,  A|-t**2*XH**2/I48.*I  XM**2U.  1**2  ) 

K(2.  2»*10.*XK«*4*6.*XJI**2 
W(2,3l*IC5.*XH**A»20.»XH**2'>3. 

«I2,  U*-C**2«XN**4/I  lt.*IXP**2«-l.)**2l 
fc«2.  J)=5.*XI«**A*12.*XH**2*15. 

NO.  2!=lC.*XP**2-*6. 

hi?.  4l=-C**2*XH»*2/I  4R.*Ixp**2U.l**21 
hi  2. 1 )  *  0. 
hi  3,  1)  *  C. 
hRITE  (6.1C1 

IC  FCRPATUX  ,25x,2E***,H-X,2tl-h  H  A  T  R  1  X,1 

1GX,3F*«*///  1 

CAU  S ESC HI  IH,2,l,C,?t,37,rrT,RANK, SOL  N 1 
AO  =  Ml.l) 

A 1  =  WI2.1) 

A2  =  VI3.11 

WRITE  I6.3C1  A,A0,Al,A2 

30  FOfcf"  AT  UX.5HA  =  , E 1 5.8 , 5X, 5HA0  =  , E 1 5 . 8, 5X ,5HA 1  =  ,El 5 . 8 , 5X , 5HA2 
1  =  «  E 1  5. 8 1 
CC  35  I  =  1.10 
CC  35  J  =  1,1] 

25  XU.Jl  =  C, 

X(l,  ll=XH**2  +  3. 

XI  1,  2M-2. 

XI 1. 21 =- 1 . 

XU,  111* AO*A 

X  I  2, 2 ) *XR**2  +  10. 

•XI 2,  21  1  • 

XI  2,  A) ■- 10. 

XI  2.  5 1  *-l . 

XI2, 51  1*A*A 1-2.*AC*A 
XI 2, 2  )  *  3.*XK**2 
XI3.21  *  3.  •»  6,  *XH**2 


X(3.6t=-t. 

XI  2,  II  »=3.«A*A2-;.*AaaC*X»**2-2.*A*A1*XM**2 
XI4,4MX»»**2*2l. 

XI  A.  *1  =  1. 

*14.7)=- 21. 

XI4.EM-1. 

XI4.  11  1=  40*4  -2.*A»A1 

X«5.4>=6.*XM**2 

XI5,  5)  =3.*Xf<**245. 

X(5.6S=+3. 

XI5.FU-*. 

XI5. 5)*-3. 

XI5.  11  )=A0*A*XM**?*A»Al*Xf'**’2-3.*A»A2 
X(fc. 5)=J.»XM**2 
X(6.< ) =15.*XH**2«3. 

X(6.«)=-2. 

Xtt. 1C >=-15. 

X(t,  11  )=AC*A*XM*»4*4.«A*A  1 *x**#4-6 .* A*A2 *XM**2 
XI 7, 7) =X0**2«36. 

X(7,E!  =♦  1  . 

XI  7.  11  1=  A*A  1 
XI fc, 7) =2 1 .*XF**2 
X(E.E)*6 .*XM**2*21. 
xte.u  =♦  6  • 

X(0,  !1  )=?.*&*A2 
XI 5, B) =10.*XP**2 
X(9, 51=13. *X6**2+1C. 

XI?.  1C  1=  ♦  1 5  . 

X(5,  11  l  =  -3.  *A1»A  *XM*»<.  +  6  .*A*A2*XM**2 
xi  ic.5  )=3.*xp**2 
X  (  1C  .  1  C  1  =  2?  .  »Xf*A*2+?  . 

XI 1C.1 n=3.*A*A2*XM»»4-2.*A*Al*XH**6 
WRITF  16.40) 

4C  FCPKATIIX  ,25X,3F»*» ,1CX, ?fhX  M  A  T  R  I  X,1 

1 CX . 3F*  »*/// ) 

CAU  SESCM1  IX..1C, 1,0,36, 37, DET.RAf.K.SCLN) 
fiO  =  X  ( 1 , 1 1 
El  =  X  12,1) 

82  =  X  I  3  ,  1 1 
B3  =  X  (4,1) 

64  =  XI5.1I 
55  =  X (6,1) 
et  =  X  (7, 1) 

B7  =  XI8.1I 
08  =  X (5,1) 

89  =  X ( 10,  1  j 

«RITE  (6,55)  BO, 01,62, 83, 84,05, B6, 67,88, B9 
55  FORMAT ( IX, 5980  =  , E l 5 . 8 , 5X , 5HB 1  =  , E 1 5 . 8 , 5X , 5HB2  =  , F 15 . 3 , 5X , 5HB3 
1=  , E  15  .8 , 5X ,  5HB4  =  ,  E  15 . 8/ / 1 X ,  5HB5  =  ,E  1  5. 8 , 5X  . 5HB6  =  ,  E 1 5 „  8 , 5X,  5H 
267  =  ,E15.8,5X,5FB8  =  , £  1 5 . 8// IX , 5HB 9  =  ,E15.8) 

CO  6C  l  =  1,15 
CO  SC  J  =  1,16 
60  YU.  J)  =  0. 

Y!  1,  11  =  120. *XM**4*144  ,*X0**2+12O. 

Y (  1,  2i=-144.*XM**2-240. 
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VI 1,21*- 240. *XM**2-144. 

VI l.*l *♦ 120. 

Yll.  :  »**72. 

Y  ( 1 «  <  1  *♦  1 20  • 

YI2..il«120.*XM**4+480.*XI'**2*e40. 
YI2.  ’1*240. *XN**2*240. 

YI2.  M»-480.*XH**2-168C. 

YI2,  :l=-240.*XM**2-480. 

Y(2,<l*-240. 

Yl 2. 71 =*640 • 

YI2, 81*4240. 

Y(2.<;1*+120. 

Yl 2,  21 =240.*XH**44?40.*X8**2 
YI3.  21  *8  40.*XM**4+48(j.*KI‘**24l20. 
YI3.‘1=-240.*XM**2 
V| 3, 51 =-480.’>XN**2-24C. 

Yl 3.4 1 =-168O.*XM**2-40C. 

Yl 3, f 1 =+120. 

YI3, 81**240. 

Yl  3 .  1C  l  =+840. 

Yl 4. 4) =120.*XM**4+10C8.*X8**2+3024. 
YI4,  51=24C.*X*I**24504. 

YI4, 41=4120. 

Y 1 4  •  7 1  =-  1008  «*XH**2— 6048 • 
Y14»£)=-240.*XM**2— ICC8. 
Yl4,S)=-240. 

Y! 4.  11  1  =  +  3G24. 

Y 1 4.  12  l=  +  504. 

Y14.  13  1  =  +  12C  • 

Y I  5 . 41 =800. *XM**4+168C . 

Yl  5.  5l=840.*XM**4  +  l6CC.*Xf'**2  4a40. 
YU,61=16H0.*XM**248CC. 

Yl  5.  71 =-l68C.*X8**2 
YU.F1  =-l60C.*XM**2-l680. 
YI5,81=-l680.*X***2-l6CC. 

YU.  10)=-168C. 

Yl  5.  121=4840. 

YU.  13  1=4800. 

YU,  141=4e4G. 

Y It,  41  =  120.*X«I**4 

YU,  51=5C4.*XH**4+240.*XM**2 

Ylb.61 =3024. *X***4+ 1008. *XK** 2+120. 

Y  V  6 , fi) =-240 . *XM**2 
YI6.91  =- 100 8 • *XM**2— 24 C  . 

Yl  6.  1C  )=-6048.*XB**2-lOCE. 

YU,  13  1  =  4120. 

Ylfc.  14  1  =  4504. 

Yl  6,  15  1=43024. 

Y I  7,  7f  = 1 20. *XM**4+ 1728. *X**» 2+7920. 
Y I  7,  8 1 =240. *XM* *2+864. 

Y(7, 51=4120. 

Yl  7,  11 1=-1728.*X8**2-1584C. 

YI7, 12 1=-240.*XM**2  -1728. 

YI7, 13 )=-240 . 

Y(£, 71=1680. *XI***4 +6048.4X8**2 
Yie.  81  =840.*XM**4+3360.*X8**2+3024. 


Y(8,5J*1680.*X***2+1680. 

Y<8.  10  )=  +  840. 

YC8,  11  )  =  -6048.  *XI<**2 
Y(8, 12 t=-3360.*XH**2-6C48. 

Y18, 13 I=-1680.*XM**2-336C. 

Y( 8. 14)=-1680. 

YIS.7>=84C.*XM**4 

Y<c,p>=1680.*XH**4*168C.«Xl<**2 

Yl5,«)=J024.*X>»**4+336C.*Xt***2+040. 

Y( 9, 1C 1=6048 „*XM**2+ 1680 • 

Y19,  12  l  =  -U>eC.*Xr**2 
Y19.  13  1  =  -3360.*XM*2-168C. 

Y19,  14  )=-6048.  *Xff**2— 3360  • 

YI9,  151---6048. 

Yl  1C  .8  1=123.  *XM**4 

Y1  10,c  1=864. *XM**4+240.*XK**2 

Y  (  10, 10  =  7920.  *X***4+1  728.  *XM**2  +  120. 

Y<  10  .1  31=-240.*X1***2 
Yl 10, ’ *1 =-1728.*XM**2-240. 

Y(  1C ,15l=-15e40.*XM**2-l 728. 

Yl  11  ,1  1)=120.*X***4+264C.*XN**2+17160. 

Yl 11 ,12)=240.*XM**2+132r. 

Yl 1 1 , 1 3 1 =+l 20. 

Yl  12,1  1)=28  80.*XA‘**4+1;8  4C.*X1'**2 
Yl  1 2,121  =  840.  *XM**4  +  5760.*Xf'**2  +  7920. 

Yl 12.1  31=1 680.  *XM**2  +  :>ebC. 

Yl 12,141 =+84 C. 

Yl 13.1 11=2024. *XM**4 

Yl 13 ,12)=3528.*XR**4+6C48.«XM**2 

Yl  13,12)  =  3C24.*XA!**4  +  7C56.*XM**2+3024. 

Yl 13 , 1*1=6048. *X8**2+3525. 

Yl 1 3 , 1  5 1  =  +  3  C  24 . 

Yl  14 ,121  =  840. *XM**4 

Yl 14 ,13l=28eC.*XH**4+168C.*XM**2 

Yl 14,14)  =7920.  *Xf'**4  +  576C.*XM**2+840, 

Yl  14 ,15)=1584C.*XM**2+288C. 

Yl  15.13)=12C.*XM**4 

Yl 15,141=1320. *X 8**4+240. *X8**2 

Yl  15  .If  »=17160.*XM**4+26  4(1.*XM**2+120. 

Yl 1, 16 )=8.*AC**2-24.*AC**2*XM**2-32. *A0*Al-4. *A1 **2~ 12. *A1* A2+4. * A 
1*B0*XR**2-2.*A*B2 

Yl  2,  16  1  =  80.  *AC*Al-48.*A0*Al*X***2+16.  *40**2*48. *AC**2*XM**2+16.*A 1 
1##2»48  .*A1*A2-2.*A*B4+4.*A*82+4.*A*B1*XM**2-4,*A*R0*XM**2 
Yl  3,  16  J  =  240.*AO*A  1*XM**2  +  16  .  *A0*  A1  *X  f<**4 +  32  ,*A0*A2-256.*A0*A2*XH»* 
12-80  .* A0**2*XM**2+ 144. *AC**2*XM**4-l  1 2 .*  Al*  A 2*1 92 .  *A  1  *A2*XM**2  -4  8. 
2»A2**2+32.*Al**2*XM**2-4.*A*85+8.*A*B2*XM**2-4.*A*BO*XM**4 
Y14, 16 )=-24.*A0**2-24.*AC**2*XM**2+96.*A0*Al*XM**2~48.*Al**2-24.*A 
1 l*42*Xf'**2-72.*AI*A2-2.*A*e7+4.*A*84-2.*A*BZ+A.*A*83*XM**2-4.*A*Bl 
2*XM**2 

Y 1 5 , 16 1=-640.*AO*A1*XM**2+256.*AO*A1*XM**4+64.*AC*A2+512.*AO*A2*XM 
1**2+ 40 C. *41*42-832.* A 1*A ?*XM**2+48 . * A0**2*XM**2- 144. *A0**2*XH**4- 1 
26 J*Al**2*XM**2+ie.*Al**2*XP**4-4.*A*B8+8.*4*B5+3.*A*B4*XM**2-8.*A* 
3B2*XR**2-4.*A*8l  *XM**4 

Yl t, 16 )=72.*A0**2*XM**4-12C.*A0**2*X«**6-384.*A0*A1*XK**4-96.*A0*A 
11*XM**6-192.*A0*A2*XM**2+768.*A0*A2*XM**4+480.*A1*A2*XN**2~552.*A1 
2*A2*XM**4-72.*Al**2*XM**4+24.*A2**2-24.*A2**2*XK**2-6.*A*B9fl2.*A* 
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-3B5*)!**«2-6.*A*82*XM**4 

VI  7,  16  )=-49  ,*A0*A1— 48  ,*A0*Al*XM**2+96,  *Al**2+4B«*Al**2*XI1**2+43  ,*A 
ll*A2+4.*A*B7-2.*A*84+4.*,S*B6*XM**2-4.*A*B3*XM**2 
VI  8.  16 1=336. *AO*A1*XM**2-272.*AO*A1*XM**4-<»6.*AO*A2-256,*AO*AZ*XM* 
i*2-l44.*Al**2*XM**2  +  112,*Al**2*Xf1**4-J36.*Al<iA2  +  lCfl8.*AJ.*A2*XM**2- 
2i44.*A2**2+8.*A*a7*X***2  +  8.*A*89-4.*A<'95-S.*A*B4*xa**2-4.*A*B3*XM* 
3*4 

VI 9, lfr)=560.*AC*Ai*XM**4-144,*A3*Al*XH**6+64.*A0*A2*XM**2-768.*A0* 
!A2*J(f'**4--1280.*Al*A2*XM**Z  +  18  72.*Al*A2*AM**4-l6.*Al**2*XM**4-96.*4 
1 l**2*Xf'**6+48«*A2**2+48.  *A2**?*XM**2  +  12.  *A*B9+12  .*A*B8*X)+**2— 12  ,*A 
3*85*XM**2-6.*A*94*XM*»4 

VI 10,lV)=17ft.*A0*Al*XM**6+8C.*A0*Al*XH**8+16C.*A0*A2*XH**4-5;2.*A0 
1*A2*XW** 6-624 •  * Al*A2*XH**4  +  5 12 .* A 1*A  2*XN**6+64.*A1**2*XH**6— 1 12  «* A 
22**2*XH**2+352.*A2**2*Xil***4+16.*A*39*X»l**2-8.*A*35*XM**4 
Vf  11  ,16)=-6C.*A1**2«24.*A1**2*XM**2-12.*A1*A2-2»*A*87-4.*A*B6*XM** 
12 

VI  12.16)=96.*Ai*‘*2*XM**2~128.<'4t»*?*XH**4+48.*Al*A2-448.*Al*A2*XM* 
1*2+48.  *A2**2-4,*A*83-8  .*  A*B7*Xrt*+2-4  •*A*B6'*XM**4 
VI  13.1£J  =  192.*Al**2*X«**4-24.*Al**2*X><**6+6  72.*Ai*»2*XK**2-132C.*A 
1 1*A2 *X 4**4-72.*A 2-**2-24.  *42**?*XM**2— 6.* A*B9— 1  2.*A+R8*XM**2-6.*A*8 
27*XM+*4 

VI  I4«lt)  =  l6  .  *A  1*42*XM**6  +  6C.*A1**2*X  8**8  +  380  .* A 1*A2*XM**4— 1024 ,  »A  1 
l*A2*Xi'**6+l6.*A2**2*X«*»2-352.=fA2»*2*XM**4-16.  *A*B9*XM**2-9.*A*8a* 
2XH++4 

Y(15.16)=256.*AI*A2*XI«**6-140.*A1M2*XM**8+88.*A2**?*XM**4-2SG.*A2 
1**?*XM**6-20.*A1*+2*XH**£-1G.*A*89*XM**4 
tePITF  (6.65) 

65  FORMAT I1X  , 29X , 3F*** , 10X , 26HY  M  A  T  o  [  X.l 

1CX,3F«**///| 

CALL  SFSCMI I Y,  15,1,0,  36, 2 7 ,0E T  RANK , SOLN ) 

CO  =  YU,  I) 

Cl  =  VI2.1) 

C2  =  V ( 3 ,1 ) 

C3  =  YI4.1) 

C4  =  Y  15,1) 

C5  =  Ylf.l) 

Ct  =  Y17.1) 

C7  =  Y  (8, 1 ) 

C8  =  Y  (9,1) 

C9  =  VIlO.ll 
C1C  =  VI  11,1) 

Cll  =  >112,1) 

C 1  2  =  Y<  13,1) 

C  1 3  =  >(  14, 1 1 
C 1 4  =  Y(  15,1) 

MFEUST/MM3/C0,CI,C2,C3,C4,C5,C6.C7,C5,CV,C19,CH,C12,C13,C1  + 

HP  IT  E  (  6  ,  NAM 3  ) 

CO  8C  I  =  1,26 
DO  8C  J  =  1,37 
6C  Z( I, J)  =  0. 

2(1,  l)=-2,*XM**2-2. 

2(1,2)  -  +  2  . 

21  l,2)=  +  2. 

2(2,2)=-2.*XM**2-12. 

ZIP, 2) =-2. 

2 1  2, 41*+ 12. 


UZ,  5>=4<, 
Z(3.2)=-<.*xi><**2 
?«  2)  =-12.  ♦XM**2-2, 
Z13,  5)=42. 

2  (  3,  <  |  =4  12. 

Z(«,4)=-2.*X^*2_30, 

Z(4,5>=-2. 

Z(4, 71=430. 

Zi4, 81=42. 

Z(  5, <)=-12. *XM**2 
Z<5, 5J=-12.*XH*#2-X2. 
Z(5,6).-12. 

Z  (  5,  6  )  =v  12. 

Z  (  5,  5 )  =4 12. 

Z  (6,5  »  =  -2.*XH**2 
z.<fc,  6)  =-3C.*XM**2~2. 

7(6, 91=42. 

Z(6, 1C  1  =  430. 

Z<7.7)=-2.,xM**2-56. 

Z  (  7 ,  8 1  *-  2  . 

Z(7,  11  )=*56. 

Z(7,  12  1=42. 

Z1 e, 7) =-30. *XM**2 

Z(8,e)--12.*XM**2-30. 

Z(8, 91  =-12.. 

Z(8,  12  1=4  30 . 

ZCfi,  13  )=«12. 

7(5,6) =-12.*XM**2 

Z(9, 5)=-30.#XH**2-i2. 
Zl9,lC)=-30. 

Z(5,  131  =  412. 

Z(9,  14  J  =  430. 

Z<  1C  ,9  t=-2.*X«*4>2 
Z(  10 ,1  Cl  =-56.*Xf<**2-2. 
2(10.141=42. 

Z*  10,1  ‘1=456. 

Z(  ll,ll)=-2. *XM**2-90 . 
2(11,121 =-2. 

Z(  11  ,1(1  =  490. 

Z> 171=42. 

4(  12,11 ) =-56.*XM**2 
Z(  12,i2)=-12«*XM**2-56. 

7(12,1 3)=-12. 

Z  (  12  , 1  71  =456. 

Z( 12,181=412. 

Z(  13.121  =  -30.*XM**2 

Z(  13,1 ■) =-3C.*XM**2-30. 
Z( 13, 14) =-30. 

Z(  13,181=430. 

Z(  13,151=430. 

Z(  14,15) =-12.*XN**2 
Z(  14,14)=-58.*xh**2-12. 

Z ( 14 ,1  * ) =-56 . 

Z(  14,151=412. 

Z(14, 20=456. 

Z( 15 , 1 4 ) =-2 . *XM**2 


2M5 

,1*1— 9C.*XM**2-2. 

£(15 

.20**2. 

2  C 15 

•  2  1 I  **90, 

7(16 

,l«ll— 2,*XK**2-132 

£(16 

,171—2. 

2(16 

,211  —  132. 

2116 

,2  2 1**2 . 

ZIU? 

i'l  <  1  —  9C  .*XM6*2 

Z(  17 

.171— 12.*XM**2-9C 

Z(  17 

.1(1—1?. 

2(17 

,2  2 1 =  +  5C . 

2(17 

,241 ■♦12. 

zne 

.171— 56.*XM**2 

Z(  18 

,lf)=-3C.*XM**2-56 

Z1  18 

.151  — JC. 

Z(  18 

,20  =  56. 

7(18 

,23)=*3C. 

2119 

,  i  S  1  —  30.*X*4**2 

2(  19 

,151— 56.*XM**2-?0 

Z(  19 

♦20—56. 

Z(  19 

,2 i I **3C • 

Z  ( 19 

.261 =  *5fc. 

Z(2C 

,10—12. *X9**2 

Z  ( 20 

,2C)=-‘>C.*XH**2-1^ 

Z(20 

.2 1 l=-90. 

Z  ( 20 

,2  6 1  =  *1 2  . 

7(20 

,2  "l  1  =+9C  . 

Z  ( 2  1 

,20—2, *XM**2 

Z  ( 2  l 

.211— l32.*XK**2-2 

7(21 

.27)  =+2. 

Z(?l 

,2ei=+n2. 

1122 

,22I=-2.*XK**2-182 

Z  (  22 

.2 ; 1 =-2  . 

Z  ( 22 

.251=182. 

Z  t  22 

• 3C 1 =+2 . 

Z(  2  3 

,22  1  —  132. *X1<**2 

2(23 

,231— 12.*XK**2-13 

Z  i  23 

, 2  4 ;  =  -1 2 . 

Z(  2  i 

,  30—122. 

Z  I  23 

.3  1 1  =  ♦  l  2  . 

Z(24 

,22»=-90.*XM**2 

Z  (  29 

,20—30. *XM**2-90 

Z(  24 

.2  5  1  =-30. 

Z  (  24 

,31) =+9C . 

7  (  24 

,  3  2  1  =  +  3  C  . 

Z  ( 25 

.241—56. *XM**2 

Z  ( 25 

.2*1— 56.*XM**?-56 

Z  l  25 

.26 1 =-56. 

Z(  25 

, 3  2 1 =  +  56 . 

Z(  25 

, 3  2 ) =  +  56 . 

Z  (  26 

.2*1 — 30.*XM**2 

Z  (  26 

,26! — 9C.*X8**2-3C 

Z  (  26 

,27! =-90. 

Z  (  26  «  3  2  1  =30  . 

Z ( 26 . ?  O  =90  . 

Z ( 27  ,261— i?.*XM**2 


.27»*-132.n>»'*2-12, 

Z<27.2f)=-132. 

ZC  27 ,34 1*12.  ■' 

Z<27,3«)»132. 

Z1 28  .27)* -2.  ♦X‘f**2 
2128.28)4-182. *XJ»**2 -2. 

Zl 28 • 3  5 ) =42 . 

Z(  28, 36)  =4-182.  '  r’  : 

2129 .25)s-2.*XH**2-20Q. 

Z129 .3C»=-2. 

Z ( 30 .25) =-l 82 • *XR**2 
Z<  30,30  =-12.  *XM**2- 182. 

2130.33)4-12. 

ZC1.3C)=-132.*XH**2 
Z  (  3 1  •  31)=-3C.*X)‘**2-132. 

2131 ,22)*-3C. 

Z132.3))=-9C.*Xf?**2 
Z ( 32 .32) =-56.4X84*2-90. 

Z ( 32  .3  2 ) =-56 . 

ZI33,32)=-56.*XM4*2 
ZI33.?3)=-9C.4XM**2-56 . 

Z  (  33  » 3 4) =  -9C . 

Z(  35.3-)=-3C.*XM**2 
Z( 35,34) =-132. 4XR4*2-3C. 

Z  1  35  » 3  5 )  =~l  32 • 

Z<  35 .?z ) =-l 2.*XM**2 
Z(  35  ,3  5J=-1B2.*X)<**2-12. 

71 35*361=— 182. 

Z( 36.33I=-2.*XM*»2 
7<  36 ,3<)=-24C.*XW4*2-2. 

Z(  1.  3?  )  =  AC*BC 

HZ,  37>=-2.4A*CC43.*AC*Bl-5.4A04B04Al*B0 
Z(  7,  -7  )=2.*A*C0*XM**24A0*82-7.*AC*B0*XM**253.*A2*B0 
Z(4. 37 J=-2.*A*Cl45.*A*C0 45.*AC*83-U .*A0*8! +7.*AC*B043.*A1*B1-5.*A 
11400 

Z(  5,  37  )=-fc.4A*C2  4f  •*A*Cl»X>'*4  24  3.*A0*84-5.*A0*82-21.*A0*Bl*XM**24l 
18.*AC*Er*XM**2-3.*Al*H2-3.*Al*b0*XM**249.*A2*8i- 15.*A2*B0»8. *AO*82 
Z  (  c  i  37  )  =  2.5A4C2*>f'**2-A.*A*CO*X.1**44AO*rt5-7.*AO'*'B2*X«*42+U.4AO*80 
14XM**5>3.*A2*B2-  13.*A2*BC*X,M**2 

Z(  7,  37  l  =  -2.*A»C3+6.*A*Cl-2.*A*C347.*A04B6-l  7.*ACi*B34n.*A0*Bl!-3.*A 
10*B045«*A1*83-11i*A1*81+7.*A1*80 

Z(  8,  27  )=-6.*A*CA4l2.*A*C2-e.AA*Cl*XM**2-2.4A*C04X)'4*24lO.*A*(;3*XM» 
1*2  46  .4An*87-3.*A04B<i-35.4AC*B3*X‘i5*2-9.*AO«B2-ll  .*A0*Br>*XM*42446.4 
2A0*B 14XK442-A14B5  4l5.»Al*e2-l 7. *A140 J.5XM442  4 10  »*A1*B0*XM**2+15.*A2 
3*8  7-33  .*A2*B  1421  .44  2*80 

Z(9. 37 )  =  -10.*A*C5 48. 4A*C24XM**249.*A*C4*XM4*2-12.*A*C1*X«**442.*A4 
1004X^4*5  43.4X0488411  „  4  A0*B  5-2  1.  *A0*B4*X')  **24  33.*  A0*B1*XM**4-13.*A0 
2*8C*XM45  42.4A0*e2*XM**2-7.*Al*B549.,*Al*B2*Xtf4*2  43.*Al*B0*XM**449. 
34A2»H5-7.4A2*R2-39.*A2*B1*XH**2438.*A24B0*XM*47 
Z(  10 ,3  71  =2.*A*C54)(M**2-4  .*A*C24XM4*4  42.*A*CO*XM44fe4AO*B94U.*AO*B2 
14xr7**5-5.*A0*B0*XM446-7.*AC*8  5*XM*  *2  43.«A2*B5-13,*A24B2*XM4*24l  7.* 
2  A2  58C*XM4*4 

7111  .37l=-2.4A*C644.4A*C3-2.4A*Cl-2  3 .*AO*B6  4  19.*AO*B3-5.*AO*B14  7.* 
1 A  1*B 6-  17.4A  1  *B 34  1 3. 4A  1*61-3.*  V*B0 
Z(12 ,37)=-6.*A*C74l2.*A4C4-16.»A4C3*XM**24Z.*A*C14XM**2-6.*A*C24l4 
1 .*A*Cfc*XM**2-25.*A0*Bl*XM*42-5.*A0*B7-49.*A0*B6*XM**2-3.*A0*B445. * 
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2A0*t)2+74.*40*a3*XK**2+Al*67+9.*Ai*84-3U*41*B3*Xr«**2-2i,  *41*32-7.* 
3Al*8C*XM»2+3fl.*Al*Bl*x***2+2l.*A2*o6-51.*A2*B3*39. *62*81-9.  *42*80 
2U3.37J»-10^*A:*Ci>20.***Ci-lC.*A*C2*XN**2*i0.*A*Ci*XK«“4H0.*A*CT 
1*X|>4*2-2C.  *A*C3*XH**4+5. 94  C*3t—3  5,  *60*87*X*t*  *2-2  5  £*40*65+5.  *40*92* 
2XH**2+55.*&0*B3*XtW>*4-35.*AC*Bl*XM?-*4+30.*A0*P4*X»»**2-5.*i:l  *83+35. 
3*Al*B5-5r*Al*fl4*XK**2  +  25.**l*81*XK*'-4-5.*Al*a0*XK**4-30.*41*B,*XK 
4**2+ 15  .**2*37-25  .*A2*84-t  5.*A2*B3*X5**2+5.*A2*B2-25.*A2*B0*X)***2  +  9 
50.*62*Sl*XH**2 

2C 14  ,371 *-l4.*A*C9+l6 ,*A*C5*X8**?-2  •  *A*C  2**X***4+6.*A»r.8*XK**2-l2 .  * 
l**C4*XP**4+6.  *A*Cl*XH**6+J5.*AO*B9-2l.*A0*88*XM**2+13.*A0*B4*XH**4 
2-5»*AO*82*XH***-15.*AO*B1*’<****6-14.*AO*3  5*XM**2-1  l. *41*89+2  l.*Al*B 
35*XH**2-9.<,Al*82*X«**4-Al*B0*XM**o+9.*A2*B8+A2*85-39.*A2*B4*:tN**2t 
45l.**2*ai*XK«*4-23.*A2*eC*****4+22.*A2*82*XM*92 
Z( 15.3  7I=2.*A*C9»XM**2-4.*A*C5*X*«*4+2.*A*C2*X<4**6-7.*A0*B9*XN**2  + 
ni.*A0*B5*AM**4-5.*AC*B2*XP**6+3.«A2*«89-l3.*62*B5*XM**2+17.*A2*B2* 
2XM**4-7.*A2*'?0*X»***‘6 

ZC  16.3 7)  =-2,*A'C  10+>».***C6-2.  *A*C 2  2  5. *A0*36-7.*AC*B3-2K*A1  *66  +  19 

1  ,*A1  »62-5.*AI«  81 

2t  17.37 l“-6.«**C 11+ 12.*A*C?-24.*A*C6*XH**2-6.*A*C4+6 ,  *A*r,i*XK**2  + 1 
»8.*A*C  10*XH**2+3  .*AO*B7+3.*A0*B4-3'’.*A0*P3*X***2  +  i;'2.*«0*B6*XK**2  + 
23.*A 1»E7-45.*A*  *86*XN:**2~ 1 5.*AI*B4+9.*A1  +B2-2 1 .*A1*B1*XH**2+66.*A1 
3*83*XH»*2-69.*.i2*86+57.*A2*B?-15.*A2*»Bl 
Z(18.37)*-10.*A*C12+20.*A*C8-8.*A*C7*Xrt**2-10.*4*C5-f..',A»C4*XM**2  + 
118  „*  A*C2*XH**4  +  14  ,*A*C  1 1  *XP** 2-26 . *A*C6*XM**4~  19«*60*6P+13.*A0*B5- 
2S.*AC*S*»*XM**2+7  7  .*AC!*86«  x!»**4-5?.  *  A0*83*XM**4*i>8.  *Al'*5? *Xf5**2  +  29. 
3*A1*  P8-19.+A 1*B7*XM**2— 45 .*Al*B5+2 1 • *A1 *62*XH**2  +4?  »*A1 *H3*XH**4— 2 
47. *A  l*E.  +  XM**4-2  .*At*84*X***2-43.*A2*B7-91.*A2*‘36*XM**2  +  23.*A2*84- 
542*6  2-  S1.*A2*81*XB**2M42  .*A2*B3*XM**2 
Z  119,37)  =-  14.*A*C13  +  2<..*A*C9  +  8.*A*C8*XM**2-18.*A»C5*XK**2+6.*A*C4* 
1XM**4+  10.*A*C12*X8**2-2C.*A*C  ”*XM*>:  4  +  10. *4*C3*Xrf**6-4l. *40*69+21.* 
2Aa*B5*X.M**2+55.  +  AC*87*XB**4-27.*A0*b4*XH**4-25.*AC*»3*XH**S>  +  14.*A0 
3*B8*XfS**2  +35  .*Al*B9  +  7 ,  *4]  *P8*X>***2+  1 3  .*A  lCii4*XH**4+l  5  .*A 
4 I*B2 *X)'+*4- 1 1  ,+A  1*B  1'  *1-70  .*41*8  55XM**!  - 1 7 .*A 2*98-65.  *A2 *87 *XH* 

5*2-1 l.*A2*B5~9.* A2*B2*XH*' 2+85. *A2*83*XM**4-S7.*A2*B1*XN**4+74.*A2 
t*64*XM**2 

Z12C.3  <)=-18,*A*C14  +  24.*»*C5*X«**2-6.*A*C5*X«**4+6.*A»C13*X‘»**2-i? 
l.*A*C8«XH**4+6.*A*C4*XK**fc  +  ?3  .*A0*83*XM**4+3 .*A0*85*XP**4- 1 5 . *A0*8 
24*XK+*Z“30. *A0*BS*XM*  +  2+  r ; .*A l*B9*XMv*2  +  Z8.*Al*B5*XM**4.+ 3 . *A1*B2*X 
3»**6+5.*A2*B9-39.*A2*B3*X4**2+51.*A2*B4*XH‘>*4-15.*A2*B2*XM**4-21 .  * 
4A2*tSl*XH**6  +  6.*A2*85*X)'**2 

1(21 ,37)=2.*A*Cl4*XM*+2-t ,*A*C ?*XM+* 4+2 . *A*C 5*X“**6+I l.*A0*B9*XM** 
14-5.  *AC*65*XH**6-13.*A2*e9*XH*-*2  +  l  7.»A2*J>*XM**4-7.*A2«B2*X4**6 
l  (22, 3  1  )  =4.  *A*C.  lC-2.*4«C6-9. *40*86+2  5.#  A  1*86-7 .  *  A1  *B3 
Z(  23,37)*  12.  *A*C  11-32  .*A»C  10*X>.**i-o-*A*C7+11.*A*C6*XB5*2  +  A0*87-5  3 
!  ,*A0*8fc*XM**2-9.*At*B7  +  7.*Al*B4-->5.-’Al*  B3'*XH**2<  94.*Al*B6*XM»*2+75 

2  .*A2*Bt-21.*A2*83 

l  ( 24  ,37)=20.*A*C12-lt.9A*Cl)tXM**2-10.*A*C8-2.*A*C7*X“**2+26.*A*C6 
1*XM*  *4-35 .*A*C10*XM**4»1 1 .*AC*B8-23 . *A0 > B7 * AH* *2 -79 . *AC*B6*XW**4-4 
23.*A  l*E8->  21  „*Al*85  +  7.*Al»e4»XB**2+69.  *A1 *06*XP**4-49 . « A1 *33 *XM**4+ 
326.*Al*B7«XM**2+4i.  .*A2*P,7-7.*A2*B4-7  t . <  A2'.'B3 * XM**2  + 1 94. *A2 *36*XH* * 
42 

i (25.371 =28. *A*C  13- 14. *A*C9- ] 4 . *A*C8*X«**2  + 14. *A*C7*XP**4-2 8.*A*C 1 
1 1*XM ** 4+ 14. *A*C4*XM**6  *  ? 1 . »A0*B9+ 7 . * A  0*6 8*XM**2-49 .*A0* B7*XM**4- 3 5 
2  ,*«C*Bi*X***6-77  ,*A  1*69  +  4?  „*A  1*8  5*XM*+2-t  35 .  *A  1  *B  7*XH**4-7  .*  A1*B4*  X 
3M**4-21.*Al*B3-vXH»*6-42.*Al*Bfi*XB**2  +  7.*A2*&B  +  !.  *A2*B5-3  5.  *A2*B4*X 
4IX**2  +  1  IS  .  *A2*66*XB**4-91  .*A 2*83*XM** 4+1 2 6. * A2*B7*XM**2 
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2  1  26  ,3  7  > «  36  .  *A-*t »4  *16  .*»*C  13*Xtt**2-2  6.**  *C9*Xfl**2  *2  .  *A*C8*X***4-2  0 
1  .***£  I  ;«XH*»4*  10  .  *A*C  7*X***6  *37.  **0*  B9«‘XM**2-1 9.  •A0*M*X****-25  .  *A 
2C*f7*X-'**6*Al*B8*Xlt**4*35.*Al*B5*A***4‘  7.»A1*B4*XR**6-U0.*A1* 
3e9*XR**Z-27.*A2*B9*7.*A2*85*XM**2*65.*A2*a7*X»t**4-A9.*A2*B4*Xit**4- 
435.*A2*83*X«**6*58.*A2*8£*»lt**2 

2  (  27 , j  7 1  *  32 .*A*C 1 4*X***2- 1 C .*  A*C9*Xft**4- 12 . * A*C1 3*XI«**4*fc . * A*C8*XH 
I **6*11  .*A0*B9*XH**4-t5.*AC*B8*XH**6-33.*AI*B9*XK**4*7.*Al*B5*XH**6 
2*51.  ■*A2*38*X***4-7.*A2*35*X**»4-21.*A2*B4*Xlt**6-10.'*A2»89*Xfe**2 
Z(2e.3  7l=-5.*AC*B9*XM**6*1/.*»?*B9*XH**4~7.*A2*65*X|***642.*A4C9*X»I 
l**6-4. *A*C14»XM**4 
2(2?  .37 ) =-9.*A 1 *B6— 2«*A*C  1C 

2 (  30.3 71—6.* ;*CU*1«.*A*C10*X8*»2»5.*A1  *87-49. *A1*86*XN**2-27.*A2 
1*86 

2(31 .371  =-lC.*A*C  12*2. *A*Cll*X***2*34.*A*C10*Xit**4*19.*Al*B8-7.*Al 
1*87*X8»»2-71.*A1*86*XM**4-13.*A2*B7-  103. ?A2«B6*XM**2 
2  (  22 .3  71  =-14.*A*C13-lC.*X*C  12*XH**2*22.*A*C1 1*XK**4*18.*A*C10-*X'4** 
16*33.*  Al*B9«-35.*Ai*88*XM**2-25.*Ai*B7*XM**4-3i ,*A1*96*XM* 

2*6*2  2*88-61  .*A2*E7*X(***2-125.  *A2*86*Xlt**4 
2(33.3  71=-18.*A*Cl4-22.*A*C13*XM»*2*10.*A*C12*XI(**44l4.*4*Cll*Xft** 
l6*77.*Al4B9*X8«*i*13.*At»e8#XH*«4-17i4Al*B7*XH**6*l5.*A2*89-l9.*A2 
2*a8*XB**2-a3.*A2*87*XK**4-49.*A2*86*X***6 
21  34,3  7>=-34.*A*Cl4*Xt***2-2.*A*C13*X***4«-10.*A*C12*XM**fc*55.*Al*B9 
I*XH**4-3.*Al*B8*XM**fc*23.«A2*B9*XM**2-4i .»A2*B8*XM**4-35 . *A2*B7*X M 
2»*fc 

Z(35.3  71=-14.*A*Ci4*XM**4*6,,*A*C13*Xf***6  +  ll.*Al*a$*XM**6+A2*B9*XM* 
1*4-2 1. »A2*B8*XM**6 

2 i 32  .37)-2.*A*C14*XH**6-7.*A2*89*XH**6 
WHITE  (t>,851 

85  ECRf'  AT  (IX  ,25X,31-***,10X, -61-2  -  ft  A  T  R  I  Xtl 

l0X.3f-***'//l 

CALL  S  E5C.H  1  (2. 36 .1.0, 36 >37. OET.  RANK  ,  SQLN  1 


DC 

7(  1 

.11 

Cl 

- 

2(2 

.1) 

02 

= 

2(3 

.11 

03 

2(4 

,1> 

04 

= 

2(5 

.1) 

05 

= 

2  ( fc 

,11 

05 

= 

2!  7 

.1) 

C7 

= 

2(8 

.1) 

08 

r 

2(5 

.1) 

09 

= 

2<  10 

,11 

01C 

= 

:*  li 

.11 

Oil 

= 

21  12 

.  1 1 

012 

= 

2  (  13 

.1) 

013 

= 

2(  14 

.11 

0)  4 

= 

2<  15 

.11 

015 

= 

2(  16 

.1) 

Dlfc 

2(  17 

,11 

Cl  7 

- 

2  (  18 

.1) 

D18 

= 

2(  19 

.11 

015 

s 

2(  20 

.1) 

020 

= 

2(21 

.1) 

021 

= 

2' 22 

,11 

022 

= 

21  23 

.1) 

02  3 

S 

2.  24 

,1) 

D24 

5 

2!  25 

,1) 

100 


025  *  2126  .11 
C26  •*  2(27  ,11 
027  >  2(28  ,1) 

028  «  2(2?  ,1) 

C29  *  Z12C  ,11 
C3C  *  2(31  ,11 
031  *  2(32  .11 
C32  »  2(33  ,11 
033  «  2(36  .11 
C34  *  2(35  .11 
035  «  2(36  ,11 

KAPf 1 1  ST /NA14/00«U1, 02,02 ,06,05,06,0  7,98 ,09, CIO, 011,012, 013, D14.D1 
15,012,017,018,  C15 ,02C,0?  1  ,C22  ,023 ,1)24 ,02 5,026 , 02  7, 02 8 ,029 , 0  3C ,  03 1  » 
2022.C32.C36.C35 
*R1T£  (5.N4K4) 

2(1.11  =  00 

2(2.  II  =  Cl  -  OC 

2(5.  II  =  C3  -  01 

2(1.21  =  C2  -  C0*XH2 

2(2.21  =  06  -  C2  -  01«X*2 

2(1.51  =  U5  -  C2»XH2 

2(7.  11  =  06  -  03 

2(5.21  *  07  -  06  -  D3*XM2 

2(3,51  =  08  -  05  -  04*XM2 

2 ( 1. 7)  =  09  -  C5*XM2 

2(9,1)  =  CIO  -  06 

2(7,  2)  =  Cll  -  07  -  Ofi*Xt>2 

2(6,5)  =  012  -  06  -  C7.XH2 

2(3.71  =  013  -  C5  -  08*X«2 


2(1.9)  = 

016  • 

■  D5*XM2 

2(11.11 

= 

015 

- 

CIO 

2 (  9.31 

= 

016 

- 

Cll 

“ 

010*XH2 

21  7.5) 

= 

017 

- 

612 

- 

C11*XM2 

2 (  5.7) 

-- 

018 

- 

C13 

- 

01 2*XM2 

2  (  3.91 

= 

019 

C  14 

m  2**82 

2 (  1.11) 

= 

02C 

- 

C14*X82 

2(  12.!) 

= 

021 

C15 

2(11.2) 

022 

- 

Clo 

- 

015*X*2 

2 (  9,5  ) 

- 

023 

- 

C17 

~ 

Clt*Xf‘2 

21  7.7) 

- 

324 

C 13 

- 

0I7*XK2 

21  5.C) 

= 

025 

CIO 

- 

D1P0XV2 

2(2.11) 

026 

C20 

0.15*X82 

21 1.  13  ) 

= 

027 

020*XM2 

2(15,1) 

02B 

C  2 1 

Z(  13.3  1 

= 

1)29 

- 

C  22 

02  1*>M2 

2(11,5) 

- 

03C 

~ 

C23 

1)2  2*  XM2 

2(5)7) 

= 

031 

“ 

024 

022*>N2 

2(7.5) 

= 

032 

025 

- 

024*X8? 

2(5.  11  ) 

- 

033 

- 

026 

D25*XR? 

2(2.  13  ) 

= 

034 

- 

C27 

026*X*2 

2(1. 15) 

- 

035 

- 

C27*Xf*2 

2(17,1)  =  -028 
2(  15  ,?  )  =  -029  -  C28*X*2 
2 (  13.5  )  =  -030  -  029*XM2 
2(  11  .7  )  =  -031  -  C30*Xf'2 
2(9,9)  =  -032  -  C31*XM2 


i 

I 


J 


I 

i 

j 

I 

! 

i 

f 


i 
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217,  in  *  -033  -  032*XN2 
Z15, 13)  *  -03%  -  033*XM2 
7(2, 151  *  -C35  -  C3**XK? 

Hi.  17)  =  -C35*XM2 
CC  ICC  L*2.18,2 
1*1-2 
XI  =  1 
LK  =  18-  t 
CO  ICO  K=  2.LK.2 
J  =  K-2 
XJ  =  J 

XJJ  =  1X.I  *  l.)/i. 

CC1=  A.C/tiXJ  ♦  i.)*(XK*»l>J+l.)1> 

CALL  !C«iT  1C. 0,1.0,. 01, XlftT.l) 

111  i  .EC.C.ANC.J.EC.2I  XIMC2  =  XINT*CC1 
IF1  I  ,f  Cw.MIC.J.fO.O)  X»KT20  =  X1MT*CC1 
Ill  I  .CC.O.ANC.J.(O.O)  XliHCC  =  XINT*CC1 
SOM  =  71  !*l,J+l)*CCi  *XIM 
C  *  C  *  SUM 
C2  =  C 

KC  CCNT  1KLC 

00=  A*(XiNTCC  -  XINT20  -  XM2*X 1NT02I 

VCSP  =  SOFT! 1  T.*A7-5.*Ar*XM2-f$G«Tl 1 5 . *A 0*XH2-3 . *A2 ) **2< 23. * AC*A2*X 
1M2  ))/(  1A.*A2*XM2  !  1 
VCRM  a  -  VU»P 

hRITF  (fc.lCCllOO.VOKP.VCPP.O? 

1CC1  FORMAT.  I  IX, SI-GO  =  #  £1  5 .8,  5X»  7HV0MP  =  ,  E 1  5.  ft  ,  SX.  7HV0RK  =  ,E1>.8,5X, 
15PC2  =  .E1S.E) 

C2CCC  =  C2/CC 

CT  =  CC*U.  -  02£00*XK*XK) 
aPITF  16, 110%)  02000,01 

1C0%  FORMAT (IX, 8H02/0C  =  , £15 . 8, 5X , l IMG < TUTAL 1  =  .F15.81 

CAU  PLOT  {  XM.XM2,  AO,  A1.A2. 80,  t)l  ,82, 33,  Ba,  85, 86,  B7,B8, 89,  CO,  Cl  ,C2, 
lC3,LA,C'),Cfa,C7,Cg,C9,CIO,Cn,Ci2,C13,Cl%,PSm,PSI21,PSI12.PSl22,P 
2SI 13,P SI  23, PSI1A.PS 12%, YC.OELY; 

GC  TC  5 
ENC 
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BIBFTC  SFSCSB  SIP.  FON.  SOLVR  OR  MATRIX  INVRTR 


SUBROUTINE  SESOM I (X, N.NB.PS ,PM ,MN2 , D,R , E) 

CCUBtE  PRECISION  NORK.SAVR ,X , Y.D.SUN, SAVEB 

CCUBLE  PRECISION  XN, XN2.AC.A1 ,42,80. B1,B2,B3 ,84,85,B6,87,»e, B9, CO , 
1C1,C2,CJ.C‘>,C5.,C£,C7,C8,C5,C1C,CII,C12,C13,C14 
DIME  NS  (CM  XI  MN1»PN2),  FORMS'))*  SAVR  (50) 

DIMENSION  (  HLO(  5C ) 

E*C  • 

R«G. 

CO  27  1*1, N 
27  SAVR ( t  )=X( 1 >  I  ) 

DO  2  1  1*1,  N 
21  IRC  (I  1*1 
IMPS)*, 4, 6 
t  NN=N+N 
NB=N 
MN=N*1 
DC  )  4  I*  1  ,N 
CO  14  J=PN,MN 

14  XU.J)=C.DO 
CO  1!  1  =  1, N 
J=t  +  N 

15  X ( I,J)=l.C0 
GC  TC  16 

A  NMN+N g 

16  JJ=NN 

SAVE  E=  X( 1 ,N  + 1 ) 

NNN=  A- 1 
C=1.CC 
CO  5  I=1,N 
KS=N-I 

IF(KK)  10,10.26 
26  LL=KK+ 1 
IJJ=1 
L=l 

WORK  =X 

DC  17  11  =  1, LI 
CO  17  J=1,LL 

IFIAeS (mCKK )-ABS(X( I  I, J)  1)18,17,17 

18  WORK  =X  ( I  I ,  J  ) 

1  =  J+ I- 1 
IJJ=J 

17  CENT  IMS 

IF(  I  JJ-1  12,2,19 

19  DC  2C  11  =  1, N 
Y=Xt  I)  ,1  ) 

X(  II  ,1  )=XI  I  I,  IJJ) 

20  X( It .IJJ  )=Y 
I Y  =  I R C (  I) 

IRD  ( I  )=  IHLDIL  ) 

IF  ID  ( L  )=  I Y 
D=-D 

2  DO  l  L  = 1 » KK 

IF(AES(X)-ABS(X(L*1,1J))7,1,1 
7  0=-0 

CO  5  J=l,JJ 


s 


V=K(  I,  J» 

Xll,J)*xtL*l,J) 

S  XtL*l,J)«V 
1  CCMUNte 
10 

IP cx  Jl  1,8,11 

11  c=c*x 

P=»H. 

- OC  12  J=l,JJ 

12  KCRK  <J  J=X(  ltJ  +  ll/x 

KK=JJ+1 
CC  3  K*1»NNN 
DO  3  J*2,KK 

3  00^5  Jji{ ^*1 ,J  *~XCK^1 • 1 ! *WORK<  J-l  ) 

5  X1M, J)=WCftK(J| 

M*=N-l 

DC  2  2  I  •  1 , N  N 
L=  I  ♦  1 

CC  2  2  J  =  L ,N 

23  «i;i!acc!,,j’,HLoijn«'»-» 

IKO(I)=lHLO(j) 

IK01J  )=  iy 
CC  25  K=1,mb 
Y=X( [,K| 

X<  1.  K)=X(J,K) 

25  X(J,K)=y 
22  f,CNT  INt6 
SUK=C,CO 
DC  28  I  =  -1 , ^ 

25  so«=suK»xu,n*s/iv!J(n 
TlrSl  =*es<  (SAvEe-SUH)/S4Vf9) 

n  KKS*-**00**”'”- 

8  F  =1. 

GO  TC  n 
CNC 


* I8FTC  I GRT 


SOBRCumE  IGRATCLL.UL.OELTA.ANS.NOEQ) 

DIMENSION  F  (  6)  ,U  (6 )  «RO) 

DATA  (MI).  1*1, tl  /  .1193C951;,-. 11930959, .33060467, -.33060469, 
A  .46623476, -.46623476  / 

CATA  CPI  II,  1*1, 3i  /  .  23355657,.  18038079. .  85662246E-1  / 

PEAL  LLIM. MULT, LL 
LLIM  =  UL 
LLIM  =  LL 
MULT  =  1.0 

IFCULIM  .C'.  LLIM  GO  TC  5 
TP  PI  *  LL  II 
LLIM  =  ULlM 
ULIM  =  TMPl 
MULT  =-1.0 
5  A  =  LLIM 

DEL  =  A8S( DELTA) 

LAST  =  1 
ANS  =  C.O 

IFCAKS  I'JLIM-LLIM)  -  .CO  CO  1) 80 , dC , L 0 
10  8  =  t  ♦  i.O*CEL 
IF(B-ULIM)40,30.2C 
20  B  =  LL  1M 
30  LAST  =  2 
4C  DO  5C  1*1,6 

X  =  (B-A)*U(  I)  ♦  .5*(A,8) 

call  meoscx.Fc  ii.noeo 

'50  CCNTINLE 

6C  ANS  =  ANS  ♦  (B-A)*(R(1  )*<F( I  )*F( 2) )  +  R I  2)*  (  F I  3 )  *F 1 4  )  I  ♦  R  C  3 ) * 
A  CFCEMFCblH 
GC  TC  (70, 30, LAST 
7C  A  =  E 
GC  TC  n 

80  ANS  =  ANS  *  MULT 
RETURN 
ENO 


AIEFTC  1M  ESC 


SUBROUTINE  I  NT  EC S < X , Y  ,N ) 
CGMMCN  Xl.XJJ 
Y  =  IX**XI)*(1.-X*X)**XJJ 
RETURN 
f  N  C 
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) 


SIBFTC  PLOTT  LIST, REF 


SUBRCL1INE  PLOT  |  XM , X«2 ,*0 • A1 , A2 , BO, SI ,B2 ,B3 ,B4, 05 »B6, B7 . B8. B9 , C 
10,C1.C2,C2.C4,C5,C6.C7,C8.C9,C10,C11  ,C12 ,C13,C14 ,SXi ,SV1  ,SX2,SY2, S 
2X3,SY3,SX4.SY4,YC,  OELY) 

DIPFNS ICN  A  A  (  1 6 )  ,  BB  ( 14 ) , CC ( 1 4 1 ,RR 1 1 3 ) ,R I C 13 > , A All  14 ) , BB1 (14 ) ,A A2 ( X 
141,052(14) 

CCUBlt  PWEC ISICN  BB.CC  ,RR  ,R  l ,  AAl  ,BB  l  ,AA2  ,BB2 
L  =  1 

XP4  =  XP2*XM2 
YF  =  YC  ♦  4.4DELY 
2  GC  TT  (5, 10, 15, 20), L 
5  CCNT  (ME 
X=SX  1 

Y  *  SY1 
GC  TC  25 

1C  CCM  INLF 
X  =  SX2 

Y  =  SY  2 
GC  TC  25 

15  CCNT  IM.E 
X  =  SX  2 

Y  =  SY  2 
GC  Tf  25 

2C  CCNT  1NLE 
X  =  SX  4 

Y  =  SY4 
23  CCNTINLE 

/AID  =  X 
A  A  ( L  +  4 )  =  Y 
L  =  L  +  1 

IF (  L.CT.4)  GO  TC  30 
GC  TC  2 
30  CCNT  [ME 
N  =  1 

22  GC  TC  (35,40)  ,N 

35  F  =  1 

2  . !  CCNT  INI.E 

Y  =  YO 

36  Y2  -  Y*Y 

GO  =  YjMA0  +  Y2*(Y2*(  (A2*XN4)*Y2t(A0*XMW.*A2*XM2))+(A2-2.*A0*XM2 

1))) 

G 1  =  Y*( ( A1-2.*AC)  +Y2*(  ( A]*X44  +  2.*A2*XM2)*Y2  ♦  ( ? . * A0*XM2-2 . *  A 1  *  X 
1P2-2,*A2) )) 

G2  =  Y  *(  ( AO-  2 ,  *  A 1  )  +  Y2*  (  2  .*A2*XM?  4-  A2)> 

G3  =  Y*(A1) 

EP(1I=  G  3 

RB 1 2  )  =  0.0 

PB  (  3  )-  G2 

PB ( 4 ) =  0.0 

EB(  5  1=  G  1 

BB(6i=  C.O 

RB(7)=  GC  -  A A ( M ) 

CO  361  I  =1,14 
361  CC( I  )  =  C.000 

CALL  RCOTSI BB,CC,6,6,RR,R  I, (ERR,AA1,BB1 , AA2 , 3B2I 

Y  =  Y  *  CECY 


> 


V. 
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IFIY.GT.YFl  GO  TC  38 
GC  TC  36 
38  H  *  F  ■*  1 

IF  (R.GT.4I  GO  TC  6C 
GO  TC  351 

40  M  =  1 
401  CO.Y  JNIF 

Y  *  YO 

41  YZ  =  Y*Y 
AF0  *  CO 

AH  «  C2  -  2.*C0*XM2 

AF2  =  C5  -  2.*C2*XM2  ♦  XM4*CC 

AF3  =  C9  -  2,*C5*XM2  +  XR4*C2 

AF4  =  C 14-  2.*C9*XM2  *  XF 4*C5 

AF5  *  -  2  .*C 1 4*XM2  +  X  F4*C  5 

AF6  =  XF4*C14 

AFB  =  Cl  -  2.*CC 

AF9  =  C4  -  2.* C2  -  2 . *C  1*XM?  +  2.*XH2*CO 

AF 10  =  C8  -  2.*C5  -  2.*C4*XF2  +  2.*XF2*C2  ♦  XX4*C1 

AF1I  =  C 13-  2.*CS  -  2.*C8*XM2  ♦  2.*XF2*C5  +  XM4*C4 

AF  12  =  -2. *04  -  2.*C13*XF2  +  2.*XH2*C9  +  XM4*C8 

AF 13  =  2 . *XM2*C 14  *  XM4*C13 

AF 14  =  C3  -  2.*C  1  +  CO 

AF  15  =  C7-2.*C4-Z.*C3*XM2+C2+2.*XM2*C1 

AF  It  =  C 12  -2.»C6-2.*C7*XF2  ♦  C5  *  2,*XM2*C4  +  XM4*C3 

AF  17  =  -2  .*03-2  .*C12*XF>2  *C9  +  2.*XF2*C8  *  XF4*C7 

AF  18  =  04  ♦  2  .*XH2*C  1 3  ♦  XF4*C.12 

AF 19  =  C6  -Z.*C3  ♦  Cl 

AF2C  =  Cll  -2.*C7  -2.*C6*XF2  «■  C4  ♦  2.*XM2*C3 
AF21  =  -2. *02  -  2.*C  1 1*XF2  +C8  +2.*XF2*C7  ♦  XM4*Cf> 

AF22  *  03  ♦  2.*XF2*C12  ♦  XF4*C11 
AF23  =  CIO  -  2.*C6  +  C 3 

AF  24  =  -2  .*C1 1  -  2.*C1C*XF2  ♦  C7  ♦  2.*XM2*C6 

AF25  =  02  ♦  2.*C11*XF2  ♦  XM4*C  10 

AF26  *  -  2.*C1C  +  Cfc 

AF27  =  Cll  ♦  2  .*XF?*C 10 

HB  (  1  1-  Y*(C1CI 

t  R  ( 2 1  =  C .  C 

FB  (  3  1=  Y *  { A  F  ? 6  ♦  Y  2*  (  AF2  7  i  ) 

88(4  l  =  C.C 

8H<5)=  Y *  1  A F 2 3  +  Y2*!AF25*\2  +  AF 24 i  ) 

8B(t  )=  C  .0 

PB(7)  =  Y*(AF19  +  Y2*(Y?*(AF22*Y2+AF21)*AF20)  ) 

88<8  !=  C.O 

C-8<9  )  =Y*(AF14  +  Y2*(Y2*IY2*(AF18*Y2  .-AF17)+AF16)+AF15)  ) 

SP(  1C)  =  0.0 

ee  ;  1  l)  =  Y*(AF8  +  Y2*(Y;.*(  Y2*(  Y2i(AF13*Y2+AFi2)+AH  1  1+AF10I  +  AF9)  ) 

DB < 1 2  I =  C.O 

B8(  1!)=Y*(AFC*Y2*| Y2*( Y2*( Y2*(Y2*( AC6*Y 2+AF 5 ! + AF4 ) ♦ AF 3 ) + AF 2 ) *AF 1 1  I 
8B  1 1  x  >  *  -AA(F  +4> 

DC  411  1  -  1,14 

411  ccm  =  o.ooo 

CALI  RCCTSI ee.CC. 13. 13.RP.FI,  IrFR,AAl,BBi,AA2,BB21 
Y=  Y  +CO.Y 
I F  ( Y  .GT.YF)  CO  TC  43 
GO  TC  41 


43  K  *  F  *  1 

IF  C  K.CT.4)  GO  TO  60 
GO  TC  40 i 
6C  N  =  X+l 

IF(N.GT.2I  GC  TO  65 
GC  TC  '2 
65  CCMINIE 
cs-l:;  =  c.19 
X  »  "0 .55 
65  CCKT  3F  liE 

T 1  =  I  1.0  ♦  XJ/t 1.0  -  XJ 

T2  *  Al/AC 

T3  =  0  .5/ 1  AC  ♦  Al) 

T  4  =  AC/Al 

IF  (T4.GE.0.0I  CC  TC  70 
T5  =  SCP  T  (-T2  1 

TFFTl  =  T3*(AL0G(T1)  +  T5*ALCGI  il.0*T5*X)/(1.0-T5*Xm 
GO  TC  75 
7C  CCNTIME 

T5  =  SCRTIT2I 

Tl-FTl  =  T  3* ( AlCG ( T 1 1  ♦  2 . C* T5*A TAN  I T 5*Xt) 

75  CC M  IME 
X2  =  X«X 

wl  =  ( 1.0-X2)*X*(X2*IX2*l86*X2  +  b3 )  ♦  Bit  ♦  BO) 

X  =  X  1  CEL  X 
CELX  =  0.20 

IF(X.G1.0.8.AND.X.lT.l.l  )  X  -  0.99 
IF(X.GT.l.O)  C-0  TC!  8C 
GC  TC  69 
PC  CCNT1M.F 
BETUFR 
ENC 


-  » I8FTC  HCOST 


SUBROUTINE  ROUTS  ( A, B, NN , NR ,PR,R I , I  ERR, A  1 ,A 2 , 6 l , B2 ) 

DIMENSION  All).  8(1),  RR (1 ) ,  RI(1),  Mill,  Bill),  A2{1>?  82(11 
DOUBLE  PRECISION  A, 8, A 1 , 8 1,A2 ,B2,RR, «I ,C,D,C 1 ,01 ,C2 ,02,F ,G , Fl.G l , F 
12,  G2  ,FFjFK«X,Y,FP 
IERR'C 
L= 1 
N»NN 
FM=1  .OC 

IF  (AN-1)  25,25.1 

1  N1=N+1 
JA=NF 

0)=CFAX1(0A8SJ A(N*1) I ,C Ae S ( B ( N*  1 )  )  ) 

C2-CAAXKCA85I  A).OABS(B)  j 
FK=N 

fH=02**i  1.00/FK)/C1**(1.CC/FK) 

IF  I  CAHSIFM)  .LT.l  .00)  GO  TC  27 
K  =  N 

CO  2  1  =  1, N 

A( I ) =A ( I )/( D1*FM**K  ) 

8 ( I ) =8 ( I ) / ( 0 1*FM**K  ) 

2  K=K-1 

AIM  )  =  A(  Ml/Cl 
8(M  )  =  E(N1)/C1 

3  CO  4  1=1, N1 
All  I  !=A(  I) 

B1U  )=  E  (  I ) 

A2U  )  =  At  I) 

4  62  *  I  )  =  F<  I) 

JA-JA-T 

NW=N-J< 

Nhl=Ah+l 
JR= J A+  1 

IF  I  JA.EC.O)  GO  TC  8 
CC  5  J=1,JA 
C--M-J 

CO  5  I  =  1  ,Nw  1 
t  =  N-J~  1*2 

aiu  i=c‘Mi  n/c 

5  tii  ( i  )=c*eu  i  )/c 

IF  I A1  If «H  )  .fC.C.OO.AMC.ai(NN+l) .E0.3.D0)  GO  TO  3 

6  X=1.C1 
Y=  .959 
11  =  1 
LL  =  1 

IF  ( Ah  .EC.l)  GO  TC  15 

7  C  =  AI 
C  =  B1 
C  1  —  A  1 

c2=-ei 

01=B  1 
C2=A  1 

00  8  I  =  1  •  NW 
F  =  X*C-Y*0+A 1 1  1  +  1 ) 

G=X*C+Y*C+B 1 ( I+l ) 

IF  (  l.EC.NH)  GO  TO  8 
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F1=X*C  1-Y*01  +  F 
F2=X *C2-Y*D2-G 
G1*a*CI+Y*CI*G 
G2=X*C24-Y*C2  +  F 
C  =  F 
C=G 
Cl  =F  1 
C2=F2 
D1=G1 
e  D2=G2 

c=c^/xi(CAes(Fi),cA8S(Gi  n 

C- CF!/>X  1(  CA8S(F2)  ,CABS(  C2J  ) 

IF  (C.LT.O)  GO  TC  10 
IF  IC.EO.LABSIFIII  GO  TO  9 
C 1-F2- FI  * ( G2/G1 ) 

IF  ICl.tC.O.CO)  CC  TO  28 
FK=(-F+F  1*(  G/Gl  I  I/Cl 
FF=-C/C1-FK*G2/G1 
GO  TC  12 

9  Cl  =G?-C-1*  <  F2/F  1 1 

IF  (Cl  .EC.O.CCI  CO  10  26 
F K  =  (  -G  +  G1«(  F/Fl  I  5/Cl 
Fh=-f /Fl-FK*f2/Fl 
GO  TC  12 

1C  IF  ; C. FC.0ABSIG2) 1  GO  TO  11 
C1=G1-C2*(F1/F2) 

IF  (Cl  .EC.O.COI  GG  TJ  28 
FH=(  -Cin‘2*(  r/F2)  I/C  1 
FK=-F/F2-FH*F1/F2 
GC  TC  12 

11  C 1  =  F 1- F2* ( G 1/G2  I 

IF  (Cl  .fC.O.CO)  CC  TC  26 
FF  =  (-MF2*(G/G2I  I/Cl 
FK=-G/02-Fh*C1/G2 

12  X=X+Fh 
Y=Y+FK 

IF  <  <X**24Y**2|.FC.<  (X+Ff-l**2  +  (Y+FKl**2t  1  GO  TO  16 

FH1=  SN  CL  ( X  +  Fh)-SXGL  1  X  I 

FK1=SNCL(Y+FKI-SFCLIY) 

IF  I FF1.NE.0..CR .FK1.NE.C.I  CC  TO  14 
GC  TC  113,141.  11 

13  LL  =  1 94 
U  =2 

14  Ll=L  L  +  1 

IF  ( IL  .GT  .2CC  I  GC  TO  16 
GO  TC  T 

15  C-A1<*2+E1**2 
X=(-A1*A1(?)-H1*R1(2)  1/0 
Y=(~A1*B’(2(+B1*A1(2' I/O 

16  IF  (JA.EC.O)  GC  TC  18 
CO  17  1*1,  N 

A2(  I  +1  )=X*A2(  I  i-Y*B2<  IH  A21  l  +  l  I 

17  82(  1  41  )*X4«?( I  I  4 Y*A2 ( 11+82(1+11 

IF  (  I  A8S(  A2IM1I  J  .GT.  1.0-E.CR  .0A8S(  B2(N+l)  I  ,GT.  1.0-81  GO  TO  21 
16  DC  2C  1J *  1 .  JR 
DC  19  1  =  1, N 


APPENDIX  B 


C  =  2.50 

XM  =  0. 40000000E  00 

A  =  0. 10000000E  01 

A0  =  0.  81325744E-02 
A2  =  -0.  972  844  16E-03 

B0  -  0. 12275631E-G2 
B2  =  -0. 72499039E-03 
B4  *=  9.53608607E-03 
B6  =  -0. 15927414E-04 
B8  =  -0. 64739868E-04 

CO  =  0. 69227686E-05 

C2  =  -0.  57445541E-05 

C4  =  0. 22511617E-06 

C6  =  -0. 13322194E-06 

C-8  =  -0. 18304866E-06 

CIO  =  0. 78165178E-08 

C12  =  0. 29131246E-07 
C14  =  0.44572204E-08 


RESULTS 

CASE  5 


A1  =  -0. 68609736E-03 

B1  =  -0. 11393347E-02 
B3  =  0.40119893E-03 
B5  =  0. 15275917E-03 
B7  =  -0. 10S73222E-03 
B9  =  -0.11241284E-04 

Cl  =  -0. 15786941E-05 
C3  =  0. 53272846E-06 
C;  0. 13238618E -05  • 

C7  =  -0. 18405079E-06 
C9  =  .  11645324E-06 

Cll  =  0. 37296582E  -07 
C13  =  C.  13214829E-07 
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DO  =  -0.78940970E-06 
D2  =  0. 12408060E-05 
D4  =  -0.36843782E-05 
D6  =  0.29335931E-05 
D8  =  0. 19504065E-05 
D10  =  -0. 11768777E-05 
D12  =  -0. 17012  411E-05 
D14  =  -0. 51517571E-07 
D16  =  0.48715251E-06 
D18  =  0. 32555807E-06 
D20  =  0. 66277360E-08 
D22  =  -0. 63921863E-07 
D24  -  -0. 83927032E-07 
D26  =  J. 67632013-08 
D28  =  0.46206709E-09 
D30  =  0. 75433856E-08 
D32  =  0. 50265173E-08 
D34  =  0. 24248057E-09 


D1  =  0. 27640558E-05 
D3  =  -0.  040000082 £-05 
D5  =  -0. 71110050E-06 
D7  =  0.40737273E-05 
D9  =  0.24398656E-06 
Dll  =  -0. 20369106-05 
D13  =  -0.  51375617E-06 
Dl-5  -  0. 25271631E-06 
D17  =  0. 61277059E-06 
D19  =  0. 93122271E-07 
D21  =  -0. 18897032E-07 
D23  =  -0. 10283680E-06 
D25  =  -0. 34505248E-07 
D27  =  -0. 49788629E-0S 
D29  =  0. 14275470E-08 
D31  =  0. 88623727E-08 
D33  -  0. 15429692E-08 
E35  =  0, 15090297E-10 


QO  =  0. 31415923E  01 
Q2  =  -0.61047082E-06 
Q2/Q0  =  -0. 19A31892E-06 


113 


VORP  =  0. 82588796E  00 
VORM  =  -0. 82583796E  00 

C  =  6.50 

XM  =  0. 15000000E  01 

A  =  0. 20000000E  01 

A0  =  0. 41087188E-02 
A2  =  -0. 19677409E-02 

B0  =  Q.38130281E-03 
B2  -  -0. 55811986E-03 
B4  -  0. 56902572E-03 
B6  =  -0.48034968E-04 
B8  =  0. 55882562E-04 

CO  =  0. 19831074E-05 
C2  =  -0.  894264  59E-06 
C4  =  0.  68938810E-06 
C6  =  -0. 17443536E-06 
C8  =  0. 39080392E-06 
CIO  =  0. 11950084E-07 
C12  =  -0.  24146492E-06 
C14  =  -0.  24739573E-06 


CASE  15 


Ai  =  -0. 19523771E-02 

B1  =  -0. 51625305E-03 
B3  =  0.27185692E-03 
35  =  0. 20099952E  -03 
B7  =  -0. 11503959E-03 
B9  =  0.21533363E-03 

Cl  =  -0.  20306899E-05 
C3  =  0.  S3055011E-06 
C5  =  -0. 20637122E-06 
C7  =  -0, 39693449E-07 
C9  =  0. 31098098E-06 
Cll  =  -0. 39021902E-07 
C13  =  -0. 40512033E-06 
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DO  =  -0. 52650429E-07 
D2  -  0. 23159584E-06 
D4  =  -0. 14998904E-05 
D6  -  0.  78764179E-06 
D8  =  0.35036516E-05 
DIO  =  -0.47970180E-06 
D12  =  -0. 53988222E-05 
D14  =  -0. 12328861E-05 
D16  =  0. 13072291E-05 
D18  =  0. 50280913E-05 ) 
D20  =  0.38015629E-06 
D22  =  -0.34301843E-06 
D24  =  -0. 25300329E-05 
D26  =  -0. 84063161E-06 
D28  =  0. 32482919E-08 
D30  =  0. 17688326E-06 
D32  -  0. 64932289E-06 
D34  =  0. 84701632E-07 

QO  =  0. 10471974E  01 
Q2  =  -0. 15016421E-07 
Q2/Q0  =  -0. 14339628E-07 


D1  =  0. 38062327E-06 
D3  =  -0. 76367925E-06 
D5  =  -0. 73885126E-08 
D7  =  0. 28239710E-05 
D9  =  0. 13379614E-05 
Dll  =  -0. 25918261E-05 
D13  =  -0. 42535094E-05 
D15  =  0. 17762290E-06 
D17  =  0. 38765723E-05 
D19  =  0. 350964  84E-05 
D21  =  -0.36900270E-07 
D23  =  -0. 13097306E-05 
D25  =  -0. 24191602E-05 
D27  =  0.  91408365E-0.7 
D29  =  0. 36867079E-07 
D31  =  0. 45445078E-06 
D33  =  0. 46538070E-06 
D35  =  -0. 52653376E-07 
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VORP  =  0. 28793975E  00 
VORM  =  -0. 28793 975E  00 

C  =  4.00 

XM  =  0. 10000000E  01 

A  =  0. 10000000E  01 

A0  =  0.  69444444E-02 
A2  =  -0. 17361111E-02 

B0  =  0.  82465277E-03 
B2  =  -0. 91145833E-03 
B4  =  0.78125000E-03 
B6  =  -0. 43402777E-04 
B8  =  -0. 13020833E-03 

CO  =  0.  53172977E-05 
C2  =  -0.  30337271E-05 
C4  =  0.14931815E-05 
C6  =  0. 14424534E-06 
C8  =  -0. 43273602E-06 
CIO  =  0. 53822889E-08 
C12  =  0. 32293733E-07 


CASE  10 


A1  =  -0. 17361111E-02 

B1  =  -0.  91145833E-03 
B3  =  0. 39062500E-03 
B4  =  0.  39062500E-03 
B7  =  -0. 13020833E-03 
B9  =  -0.43402777E-04 

Cl  =  -0. 30337271E-05 
C3  =  0.  86769225E-06 
C5  =  0.  86769225E-06 
C7  =  -0. 43273602E-06 
C9  =  -0. 14424534E-06 
Cll  =  0, 21529155E-07 
C13  =  0.21529155E-C7 


C15  =  0. 53822889E-08 


>' 
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DO  = 

-0. 27081060E-06 

D1  =  0. 14618229E-05 

D2  = 

0. 85993357E-06 

D3  =  -0. 24770679E-05 

D4  = 

-0. 38892866E-05 

D5  =  -0. 14446701E-05 

D6  = 

0. 21444442E-05 

D7  =  0. 55700687E-05 

D8  = 

0. 51091772E-05 

D9  =  0. 13936700E-05 

DIO 

=  “0. 10670710E-05 

Dll  =  -0„  38699997E-08 

D12 

=  -0. 57397364E-05 

D13  =  -0. 30777047E-03 

D14 

=  -0. 78396422E-06 

D15  =  0. 31178924E-06 

D16 

=  0. 14493480E-05 

D17  =  0.  28736014E-05 

D18 

=  3. 23597631E-05 

D19  =  0. 16021832E-05 

D20 

=  0. 22322972E-06 

D21  =  -0. 4668c O16E-07 

D22 

=  26815118E-06 

D23  =  -0. 64046568E-06 

D24 

=  -0. 81407119E-06 

D25  =  -0. 58064111E-06 

D26 

=  -0.22029153E-06 

D27  =  -0.  34721103E-07 

D28 

=  0. 26964840E-08 

D29  =  0. 18191922E-07 

D30 

=  0. 52525373E-07 

D31  0. 84124962E-07 

D32 

=  0. 80707635E-07 

L»33  =  0.462“  185E-07 

D34 

=  0. 14774596E-07 

D35  =  0.20130^  J7E-08 

QO  =  0. 15707962E  01 
Q2  =  -0. 12776789E-06 
Qg/QO  =  -0.  81339574E-07 
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VQRP  =  0.42923705E  00 
VORM  =  -0.42923705E  00 


APPENDIX  C 


LIST  OF  SYMBOLS 


A  =  constant  = 


Aj  =  constants 
bj  =  constants 


C  =  constant  = 


2 (nr  +  1) 


i  =  0,1,2 
i  =  0,1,2,.. .,9 


Ga 

MW0 


C.  =  constants 

l 


d.  =  constants 


i  =  0,1,2, ,  14 
i  =  0,1,2, ... ,35 
e  =  eccentricity  of  the  cross-sectional  ellipse 
et  =  constant 

=  components  of  the  body  forces 
G  =  constant  =  mean  pressure  gradient 


K  =  constant  = 


2W§a3 

Rv 


m  =  constant  related  to  the  eccentricity  by  e 


_ 1  n 


m 


m  1  and  e  =  —  \Tl"  -  m2  when  0  <  m  s  1 
m 

P  =  pressure 

,  <jy  =  dimensions li zed  velocity  components 

=  dimensionalized  velocity  in  the  x-direction 
r,0,y  -  cylindrical  coordinates 

R  =  radius  of  curvature  of  the  center  of  the  pipe 
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W  =  nondimensionalized  velocity  in  the  9-direction 
W0  =  velocity  along  the  central  axis  for  the  flow  through  a  straight 
elliptic  pipe 

Wt  =  first-order  approximation  of  W 
VV2  =  second-order  approximation  of  W 
H  ~  absolute  viscosity 

u 

v  -  —  =  kinematic  viscosity 
P 

p  =  fluid  density 
ip  =  stream  function 

ibt  =  first-order  approximation  of  the  stream  function 
4>2  -  second-order  approximation  of  the  ctream  function 
The  '  superscripts  with  parameters  denote  dimensionalized  variables. 
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